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Section  I 


INTRODUCTION 


An  objective  of  the  present  research  program  was  to  develop  a  finite-element  based 
procedure  for  analysis  of  free-edge  delamination  specimens  using  through-the-thickness 
elements  and  including  l^oih  stretching  and  trending  effects.  This  necessitates  the  use  of 
muliilaver  plate  themies  which  can  simultaneously  consider  bending  and  stretching. 
Multilayer  plate  theories  have  l)een  developed  using  assumptions  on  displacements  or  on 
stresses.  I  he  former  class  of  theories  may  be  classified  into  two  groups,  vb-, 

1.  Theories  based  upon  assumed  variation  of  displacement  as  a  polynomial  in  the 
transverse  coordinate  over  the  entire  thickness. 

2.  Theories  based  on  assumption  of  piecewise  linear  variation  of  displacements 
over  the  thickness  with  'nodes'  at  the  interfacial  surfaces. 

The  first  group  of  the  displacement  theories  has  been  found  to  be  inadequate  for  repre- 
senution  of  behavior  of  composite  laminates  where  the  material  properties  along  fiber 
are  significantly  different  from  those  in  the  directions  across  fiber.  The  second  group 
of  theories,  culled  the  discrete  laminate  theories,  apparently  is  a  better  candidate  for 
further  consideration.  These  discrete  laminate  theories  have  been  described  by  Srinivas 
[l]  and  Sun  and  Whitney  [2],  among  others,  and  solution  schemes  have  been  proposed. 
However,  for  arbitrary  geometry  and  a  larg''  number  of  layers  one  has  to  resort  to 
numerical  procedures.  The  finite  element  method  has  been  extensively  used  for  the 
analysis  of  plates  (  e.g.  Reddy  [3]  and  [4],  Davis  [5]  among  others). 

Davis  and  Mawenya  [5]  proposed  a  general  finite  element  formulation  using 
quadratic,  isoparametric,  multilayer  plate  elements  which  allowed  layers  to  deform 
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locally  with  no  restriction  imposed  on  the  relative  properties  of  the  constituent  layers. 
Since  this  formulation  considers  the  transverse  shear  deformation  in  all  layers,  it  is 
applicable  to  any  arbitrarily  layered  plate.  However,  the  stresses  are  discontinuous  across 
the  interfaces. 

As  a  starting  point  in  the  present  research  program,  Davis  and  Mawenya's  approach 
was  used  to  develop  a  finite  element  solution  for  analysis  of  laminated  plates.  Davis 
and  Mawenya  did  not  give  details  of  the  theoretical  and  numerical  formulation  they 
used.  In  this  report,  the  theory  is  re.stated  in  variational  form,  and  an  extension  oi 
the  general  variational  theory  is  specialized  for  implementation  in  u  finite  element 

computer  program.  Its  effectiveness  for  analysis  of  plates  and  also  its  inadequacy  in 

modelling  stresses  in  free-edge  delamination  specimens  are  noted. 

Section  II  contains  a  summary  of  the  equations  governing  bending  and  stretching 
of  laminated  plates.  In  this  section,  the  kinematic,  equilibrium,  and  constitutive 
equations  for  a  lamina  are  derived  based  on  an  extension  of  Mindlin's  theory  of  plates. 
The  displacement  field  is  assumed  such  that  the  rotation  of  each  lamina  is  an 

additional  variable  independent  of  transverse  deflection.  The  set  of  coupled  field 
equations  and  interlaminar  continuity  conditions  is  written  as  a  self-adjoint  matrix  of 
operators.  Consistent  boundary  conditions  are  identified,  and  a  general  variational 
formulation  for  the  purpose  of  finite  element  approximation  to  the  problem  is 
developed.  Section  Ill  discusses  the  finite  element  formulation  and  computer 
implementation  of  the  theory.  Some  illustrative  examples  and  comparisons  of  results 
against  some  alternative  solution  schemes  are  discussed  in  section  IV.  This  section  also 
includes  application  of  the  procedure  developed  to  analyse  a  multi-layer  free-edge 

delamination  specimen. 
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Section  II 


BENDING  AND  STRETCHING  OF  LAMINATED  PLATES 

2.1  INTRODUCTION 

Structural  elements  composed  of  an  arbitrary  number  of  orthotropic  layers  can  be 
approximated  by  finite  element  pnKedures.  In  these  composite  elements,  each  layer  may 
have  different  thickness  and/or  elastic  properties  and  different  orientation  of  axes  of 
material  symmetry.  In  the  following,  we  summarize  the  governing  equations  for 
bending  of  plates  based  on  the  following  assumptions. 

1.  Loads  are  carried  primarily  by  bending  and  stretching  of  the  plate. 

2.  Sliding  of  one  layer  past  another  is  impossible. 

3.  Plane  sections  normal  to  the  undeformed  surface  of  each  layer  remain  plane 
but  not  necessarily  normal  in  the  deformed  configuration  and  the  in-plane  dis¬ 
placements  vary  linearly  over  each  layer. 

4.  The  transverse  displacement  is  independent  of  the  transverse  coordinate  i.e.  the 
transverse  strain  vanishes. 

5.  Deformations  and  rotations  are  small  and  the  material  is  linear  elastic  so  that 
the  linear  theory  of  elasticity  is  applicable. 

Since  the  neutral  axis  is  a  priori  unknown,  bending  and  stretching  are  coupled  with 
respect  to  an  arbitrary  plane  of  reference. 

With  the  above  assumptions,  using  the  rectangular  Cartesian  reference  frame,  the 
kinematic  field  variables  consist  of  three  displacement  components  for  an  arbitrary 
px)int  on  a  reference  surface  defined  by  a  constant  value  of  the  transverse  coordinate 
in  the  reference  (undeformed)  configuration  along  with  the  values  of  the  rotations  of 


segments  of  the  ray  along  the  transverse  axis  througli  this  point  defined  by  intersection 
with  interlaminar  surfaces.  The  numixr  of  field  variables,  therefore,  is  2m+3  where 
m  is  the  total  number  of  layers. 

2.2  EQUATIONS  GOVERNING  BENDING  AND  STRETCHING  OF  LAMINATED 
PLATES 

2.2.1  Introduction 

The  generalized  equilibrium  equations  represent  the  integral,  over  the  thickness,  of 
the  three-dimensional  equilibrium  equations  and  of  the  first  moment  of  the  equilibrium 
equation.  The  constitutive  equations  are  stated  for  a  linear  elastic  monoclinic  material. 
I'or  implementation  in  a  Ritz  type  finite  element  appro.ximation  procedure,  the  problem 
is  formulated  as  a  set  of  self-adjoint  field  equations  with  consistent  boundary  operators. 
The  index  notation  is  used  throughout.  Latin  indices  uke  on  the  range  of  values  1,  2, 
and  3  whereas  Greek  indices  take  values  1  and  2.  Subscripts  following  a  subscripted 
comma  denote  partial  differentiation  with  respect  to  the  coordinates  defined  by  the 
subscripts.  Summation  on  repeated  indices  is  implied  except  where  indicated  otherwise. 
A  pair  of  indices  within  parentheses  denotes  the  symmetric  part  of  the  tensor  described 
by  the  .subscripts  and  a  .single  super-  or  subscript  within  parentheses  denotes  'no  sum' 
on  that  index. 

The  actual  displacement  vector  at  any  point  is  a  function  of  the  coordinates  (jc,) 
of  the  plate.  Assumption  of  transverse  displacement  being  independent  of  the 
transverse  coordinate  x,  makes  it  a  function  of  (x„)  only. 
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2.2.2  Kinematics 


('omponenls  of  llie  displacement  vector  lor  each  lamina,  assuming  linear  variation 
ol  inplane  displacement  over  the  thickness  of  the  lamina  (Fig.  l),  can  be  stated  in  the 


a  I  a 


u^\x.)  =  u!^\x  )  =  uAx  ) 

Here  a  rectangular  ('artesian  frame  of  reference  is  used.  are  components  of  the 

displacement  vector,  and  v*'*  are  the  'inplane'  displacements  at  the  kth  interface. 
are  the  components  of  the  rotation  of  the  kth  layer  in  the  a  — 3  planes.  For  infinites- 

Qu 

imal  elastic  deformation  i.e.,  — -<<1,  the  strain-displacement  relationship  is: 

QXj 

€  ,  =  4-(“  )  =  U,  (2) 


— V**  .V 

•j  2  *'•' 


Theref  ore, 


where 


M)  _  U)  ,  .  (<)  U) 

e"  =  1(„“'  +«<") 
a3  2 


4’  =  =  0 


e“'  =  v“’  =i(v'‘’ +v'‘>) 


U)  ,U)  _  1  /  ,U)  ,  .(O-v 
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2.2.3  Equilibrium  Equations 


The  three-dimensional  equations  of  motion  including  body  forces  are: 

U)  ,  Ai)  U)  -it) 

(T  =p  u.  15; 

where  cr,^  are  components  of  the  Cauchy  stress  tensor  and  /,  those  of  the  body  force 
per  unit  volume.  Equation  (5)  is  difficult  to  satisfy  exactly.  To  eliminate  the 
dependence  on  the  coordinate  Xj,  (5)  can  be  restated  in  the  form 


+  u''yx‘‘^y‘dx^' =  0,  n=0.1,...oo.  (6) 

(i 

As  an  approximation,  in  general,  (6)  :s  enforced  only  for  n=0,l.  Higher  order  theories 
would  use  higher  values  on  n  as  \xell.  Equation(6)  for  n=0  is: 


a.  I'or  i=l,2,  setting  +  following  (l): 


.,«)  *(i)  .  r.(^>  „U)  ••(/I) 

N..+  ia’-,-<r^)  +  F  -F  v  —R  (f>  =0 

o0,/3  o3  a3  o  « 


where 


k 

.,(<)  r  o;  ,  (n 

•?) 


'I 

(f  ,/?  )  =  Ip  (1  ,x,  )<fx^ 


(10) 


b.  For  i=3. 


k 

//  U)  ,  «)  .  A  (<)  ••(O'l  J  O)  ,, 

^«^o3.o+^33..-.+^-P  ^^^^3  =‘* 


+(o-^*  -  o-:;''  ’)+  /■‘f  -  -0 

*’of,or  33  33  3 


where 


gI!^=  f 

*f) 

,.(n  r  Ji( )  j  (< I 

^'  l  =  J  /3  ^^3 


For  n=l,  i.e^  taking  the  first  moment,  considering  i=l,2. 


KAfl+‘^«3.3+4-#‘  ''a  -P  (^3  ^‘^*3 


j.*  At),^U)  r/k) -ii)  ,<0  va)_^ 

^  <^a3  “  Q.  ^  ''a  “ 


where 


■1 


(O  (X)j  (i) 

O'o^^3  <^^3 


.U).Ji) 

3  ^^*3 


]'p'‘>(x';Vdx‘; 


In  the  approximate  the.  ry  considered,  the  first  moment  equation  for  i=3  is  ignored. 
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To  eliminate  the  time  derivatives  in  (7),  (11),  and  (14)  and  to  include  the  initial 
conditions,  convolution  of  both  sides  of  the  equations  with,  g(t)=t,  a  function  of  time 
is  performed  [8]  and  [9].  Kquations  (7),  (ll),  and  (14)  can  then  be  written  in  the 
following  form; 

g*  + t^g*  >  -  g»  d:  '+g*  g;:  v  v-;  ’  -  /“  v  = o  (20) 

where  (*)  represents  the  convolution  prcxluct.  l-'valuating  the  convolution  integrals  for 
terms  involving  time  derivatives 


'(///Uj-) 

g*  J  gil-r)  iV^\x^,T)dT 


”U), 

w  dr 


=  -  -  >v‘'\xp,0) 


where  »v‘*^x^O)  and  w‘*\xg,0)  are,  respectively,  the  initial  conditions  for  transverse 
velocity  and  displacement.  Similarly, 


^(k)  ^  .pU)[  _  ,H,<*>(jc^,o)-w‘'  \x^,0)] 

y{k)  ^  \xp,0)]  --  - 1 


Substituting  (24)  through  (26)  into  (IS)  through  (20): 


d 


(27) 


-UK,  ,.«)  u).  „!/'  ,, 

-a33  )  +  g*/', W  +X  =0 


-1_  _*  L-'"  o'".,'"  o'"j.'"  _L  v'‘‘^  _  O 

*  ,-(X)  ,.  _:t*U)  ^ylU^.,^^U)  „«)(/) 

+  ^ 


.(/)  ,yUlU)  ,yU),U)  ,  „(0 


(28) 

(29) 


Eqs.  (27)  through  (29)  are  the  spatial  equilibrium  equations  for  the  motion  of  each 
lamina  of  the  laminated  plate  in  terms  of  laminar  force  resultants. 

2.2.4  Constitutive  Equations 

2.2.4.1  Stress-Strain  Relations  for  Linear  Elastic  Materials 

A  material  is  said  to  be  ideally  elastic  if  the  material  completely  recovers  its 
original  shape  upon  the  removal  of  the  forces  causing  the  deformation.  7'he  generalized 
Hooke's  law  relates  the  nine  components  of  stress  and  the  nine  components  of  strain 
by  a  linear  relation.  Assuming  an  initially  unstressed  reference  configuration,  and  a 
rectangular  Cartesian  reference  frame,  this  can  be  expressed  as: 


O’.  ,  = 

ij  ijkl  kl 


(30) 


where  are  components  of  fourth  rank  isothermal  elasticity  tensor.  Owing  to 

symmetry  of  or,^,  i.e.,  in  the  absence  of  body  couples,  we  have: 


^ijkt  ^  jiU 


Furthermore  since  €*,  =  €,^: 


E.  t,  -  E. ... 

tjk!  ij/t 


If  a  strain  energy  function  exists  then 


^ijki  ^tiij 


(31) 


(32) 


(33) 


For  a  monoclinic  system  with  two  ortht^onal  planes  of  symmetry  the  stress-strain 
relations  for  any  layer  k  in  the  reduced  form  with  respect  to  a  global  plane  of 
reference  (second  rank  tensors  written  in  vector  form)  are: 
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,.(X) 

^1111 

jjU) 

'"1122 

Ak) 

^1133 

0 

0 

XX) 

1112 

11 

^^22 

.,(X) 

'"2211 

,.(X) 

"2222 

p(X ) 

^2233 

0 

0 

,Xx) 

2212 

(X) 

*"22 

^U) 

jpik) 

0 

0 

XX) 

(X) 

33 

3311 

3322 

3333 

^3312 

33 

<^23 

0 

0 

0 

Ak) 

^2323 

Ak) 

^2313 

0 

2£‘'^ 

"23 

<^13 

0 

0 

0 

Ak) 

^1323 

^k) 

^1313 

0 

13 

*^12 

„(X) 

•^1211 

^k) 

'^1222 

^X) 

■^1233 

0 

0 

Ak) 

^1212 

2£‘^’ 
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Ikluation  (34)  can  be  written  in  the  indicia!  form  as: 


U)  _  .,U>  U)  Ut)  (X) 


(X)  _  (X  )  (X) 

tr  ,  =  2/i  ,  ,e  , 

aj  a3y3  >3 


(/)  _  ,.U)  U)  ,  r?(^)  U) 
33  33>o  yl  3333  33 


Solving  (37)  for  €33^  and  substituting  into  (35) 


(34) 


(35) 

(36) 

(37) 


U)  _ 

01  afiyh^yb 


(38) 


where 


Ut) 

■j^t)  (i)  ^01033  p(*) 

a0yb~  EaPy()~  „(X)  '®33y6 

■^3333 


(39) 


and 


p(X) 

■WX)  _  ^a<333 
^0033  Ut) 

^3333 


(40) 


Substituting  (3)  into  (38),  (36),  and  (37>. 


■W*)  (k)  (k)  (X) 

afiy^yt  afiyb  ^3  ^ofi 


ik) 


^0033  ^33 


(41) 


(42) 


.(X)  _  ^x) 


JX)  ^  ^X)  Jx)  Jk)  ,  ^x) 


JX) 


*^33  ^33>8S6  ^33y6'*^3  ^3333^33 


(43) 


We  note  here  that  the  assumption  of  ix,*’  constant  over  the  thickness  of  each  lamina. 


II 


implies  u“ ’(x',)  =  €‘,V  =  ().  This  wouUl  result  in  the  terms  containing  e\\' 

dropping  from  (43).  We  also  note  that  this  would  mean  that  (37)  and  (35)  would 
not  contain  €jj’  and,  therefore,  there  would  be  no  question  of  eliminating  633^  between 
these  two  equations.  However,  most  plate  theories  suffer  from  this  defect  viz.  e,*’  is 
eliminated  in  the  constitutive  equations  (35)  through  (37)  but  is  set  equal  to  zero  in 
the  equation  for  (Tj,  i.e.,  (43). 
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(50) 


A/  =  7  A/'^; 

(\fi 

k  1 


(51) 


<•1 

w  here  T,*,  is  the  distance  of  the  center  of  the  kth  layer  from  the  point  of  application 
of  the  resultant 


2.2.5  The  Interlaminar  Continuity  Equations 

for  tlie  continuitN  oi  ir.ictions  and  displ.icemcnis  to  be  satisfied 


U  )  .  If,  Ifll 

V  +  /,  0  =  V 


and 


(-0  _ 

33  " 


<r 


33 


where 


(52) 

(53) 

(54) 


(55) 


=  components  of  the  shearing  stress  at  the  top  of  the  kth  layer. 

=  components  of  the  shearing  stress  at  the  bottom  of  the  (k+lXh  layer, 
v'*’  =  components  of  the  m-plane  displacement  for  the  kth  interface. 

0l*’  =  components  of  the  rotation  of  the  kth  layer. 


2.2.6  Summary  of  Field  Equations  in  Convolution  Form 

For  self-adjointness  of  the  set  of  operators  consisting  of  those  appearing  in  the 
kinematics  equations,  the  constitutive  relations,  the  equilibrium  equations,  and  the 
continuity  equations,  noting  that  the  equilibrium  equations  had  to  be  transformed  (27) 
through  (29),  it  is  necessarv  to  espress  the  remaining  field  equations  in  convolution 


1.1 


form  as  well.  i.e.. 


a. 

1-quilibrium  Ikiuations  ((27)  through  (29)); 

s  Q...  +  «  -‘^13  )  +  g'  vv-  +  JV  =0 

(56) 

*  _L  ^9)x,  »  ,  .(/C)  „«)  U)  m(<  )  j  <' )  .  v*' >  /^ 

(57) 

+  ''o  “■'  <^0  =0 

(58) 

b. 

Kinematic  Relations  ((2)  through  (4)) 

(59) 

s*C  = 

(60) 

vOg‘(vv;^’+<') 

(61  ) 

c. 

(Constitutive  I'-quations  ((44),  (45)  and  (47)) 

^^0,^  “  ^\py(,eyb'^  r^aftybKyb 

(62) 

f  ^  D^')  (^')  i  9') 

«  ~  o  ^ofiybfyb'^  ^^afiybKyb 

(63) 

(64) 

d. 

Continuity  Equations  ((52)  through  (55)) 

*_♦<*) _ u*i) 

(65) 

(66) 

«)  *  «'I) 

g*w  =  g*w 

(67) 

«*‘^33  =S^‘^33 

(68) 

The 

field  operators,  for  layer  k,  in  self-adjoint  form  are: 
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-/■>“’  0 


0 

-g*L  g* 
0 


[Af '  = 


»L  -R 


()  0 


*  .(<)  (n 

-S  \pyi  g*  0  -g*/i..p>6  O  0 


0 

0 


e  o/3y6  0 


0 

0 

0 

0 


0 


0 


0  {) 
g*L  0 


0 

0 

0 

0 


0 

-g*i^  g* 
0  0 

0  0 

-g*  0 


U) 


0 

0 

0 

0 


0 

-F 

0 

-f* 


,(i) 


Qa 


0 

0 
0 
0 

0 
0 

0 
-g*G' 
g 


it) 


*  0 


0 

0 

0 

-g* 

0 

0 

*  da 
g* 


(69) 


where 


L  =  1(8  -4-+S  .-^) 

2  V"  a/3  Qa 


The  domain  of  the  operators  is  the  direct-sum  space,  under  the  convolution  product,  of 
the  spaces  consisting  of  admissible  ,Q^^\  in  that 

order.  The  operators  in  (69)  are  the  ones  in  (57),  (62),  (59),  (58),  (63),  (60),  (56),  (64), 
and  (61)  excluding  those  associated  with  interlaminar  tractions.  The  displacement 
continuity  (66)  and  (67)  can  be  directly  incorporated  in  the  system  of  equations.  At 
the  same  time,  addition  of  operators  associated  with  the  interlaminar  traction  completes 
the  set  of  field  equations  for  the  system.  Traction  continuity  ((65)  and  (68))  is 
implicitly  satisfied  by  using  only  one  set  of  tractions  for  each  interface.  I'he  system 
set  of  field  equations  then  has  the  form: 


[0] 

[d^Y 

[0] 

(0} 

^pf2) 

{0} 

lO] 

'Y  (01 

*>} 

10} 

[c"’] 

ird'+W-') 

(70) 
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or,  symbolically. 


[X]{y)  = 

Here 


[5^‘Y  = 

g*  0  0  0  0  0  0  0 

0  000  00g»00 

= 

V 000000  00 

0  00000  -g»  00 

(71) 

(72) 


I  (ihT  .  (/)  (/)  .,u)  ,0)  u<  .jt)  (li  ^  (/ ■ 

iu  }  =  <v  ,e  ,K  ,2€  ,.G  > 

a  op  op  (kp  <»p  o3 


(73) 


{a-  }  =  <ct_3  ,a,3  > 


(74) 


o  O 


.(/() 


.O.O.-g’/','- A"'", (),()> 


(75) 


{()!'  =  <0,0  >  (7t,) 

{cr^Y  =  <g*a^3\0,0,0,0,0,g*cr3*3‘\0.0>  (77) 

{£rY=  <-^cr*^3^\0.0,-/^gV3^\0.0,-^<733'’’.D,0>  (78) 


[0]  = 


0  0 
0  0 


(79) 


Continuity  of  tractions  is  ensured  explicitly  by  using  the  interfacia]  traction  as  the 
field  variable  [o’'**]  in  the  manner  expressed  by  (70),  i.e.,  o-',*'  =  0",^  and  do  not 

appear  as  field  variables.  The  operators  and  their  adjoinls  ,[C‘*T 

represent  coupling  between  field  equations  for  the  layers  and  the  continuity  of 
interlayer  displacements.  Explicitly,  for  the  interface  between  the  kth  and  the  (k+l)th 
layer,  these  have  the  form 


lb 


1.4/*’ 


«-  0 
0  0 
0  0 
»^*0 
0  0 


0  0  I^g-  0 

0  0  0  0  o 


lol 


0 

0 

0 

0 

0  0  0  0  0 

0  0  0  0 
r' 
o 
0 
0 
0 
0 
0 
0 
0 
0 
0 


lol 


0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

s' 

0 

0 

0 

0 


0  0  0  0 
0  0  0  0  0 


O  0 

o  >. 


I.AI**-' 


0  0 
0  0 


(il 

V 

U) 

y*  *  ^ 
'ufl 

.(*) 

K 

^U) 

|m' 


.<*) 

(k) 


33 

(itw) 

y 

(X  tl) 


.,(X4l) 


,(X4l) 


(X*l) 


Ia; 


(X4l) 


ad 

(X4l) 


|2€' 


.(i:4l)i 

03 


,(*♦1) 

o3 

,(*+!) 

o3 


.  ..(*•)  ..(S) 
K 


-g*C^*’-Z**’ 


0 
0 
0 
0 

^  ,^S.I)  J  (i-.Di 

|0 

|o 

lo 

|o 

0 
0 
0 
0 


2.2.7  Self-adjointness  of  the  Operator  Matrix 

For  the  operator  matrix  [X]  to  be  self-adjoint,  a  sufficient  condition  is  that  the 
elements  of  X  satisfy  the  relationship 

<  /, .  8 j  >  A-  =  <g, .  Xj,  /,  >  A  Boundary  terms 

where  <  ,  >,  is  a  bilinear  mapping  over  the  space  of  functions  defined  over  the 
region  R.  The  elements  of  X  satisfy  this  requirement  in  the  sense  of  inner  product 
defined  as 

<f,g>  =  f  f(x)g(,x)dR 

Specifically,  if  the  operator  matrix  /4'*’  is  self-adjoint  and  ([^“’],  [5‘*’]0:  ([C**!  [C/*’]^) 
constitute  adjoint  pairs,  (X’)  is  self-adjoint.  Considering  A,j  and  A^,  for  the  kth  layer, 
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based  on  Green’s  theorem  (Kreyszig  [1979]); 


<  F* >  =— e*v*^  > 


afi'o  ’'(a.B)  ,^1=^ 

M)_  *.»), 


(O  *»,U). 


+  <Kfi‘n,^g*K>,uy+  <<^g*K;'np  >.*, 


Similarly  for  and 


<0„  'i*  ^  ap,p^  jjny 


,(/) 


(/)  .  ,jn 


.1^  .». 


xU) 


and  for  and  Af/: 


<  ^''\g*QZ  >y*)  ="  <c« 


+  <  > 


4*’ 


..W) 


(80) 


(81) 


(82) 


The  remaining  elements  of  constitute  algebraic  operators  which  are  self-adjoint  or 
consist  of  adjoint  pairs.  The  adjointness  of  and  [//*T.  and  of  [C‘*’l  and  [C**T  is 
obvious. 


2^8  Consistent  Boundary  Conditions 

Referring  to  Sandhu  [1975],  consistent  boundary  conditions  for  the  problem  are: 


on 

on 

on 

4” 

(83) 

on 

5s’ 
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where  and  C*/'  are  the  consistent  boundary  operators  and 

1  2  3  4  5  6 

^  =  <f> 


lixplicitly, 

^,U) 

„(/) 

=  ^.s 

^'2 

II 

o 

- 

(Corresponding  internal  jump  discontinuity  conditions  are: 


=- 

■g*(gV) 

on 

M) 

^1, 

g*id^^M%y  =- 

on 

M) 

^3, 

_  _ 

on 

on 

M) 

^2i 

g*(gV) 

on 

sV 

At 

II 

on 

sV 

(84) 


Here  ...  ,S^^'  represent  surfaces  imbedded  in  the  interior  of  the  region  R,  A 

prime  over  any  quantity  denotes  the  jump  in  that  quantity  across  the  surface  of 
interest,  e.g. 


/  =  f*-r 


where  ±  denote  the  two  sides  of  the  interior  surface.  Quantities  i=l,  2, 


6, 


denote  the  specified  values  of  the  jump  discontinuities. 


23  VARIATIONAL  FORMULATION  OF  THE  PROBLEM 


23.1  The  Basic  Variational  Formulation 

Based  on  (A.20),  the  self-adjoint  form  of  the  field  equations  given  in  (70)  the 
boundary  conditions  (83),  and  the  jump  discontinuity  conditions  (84),  the  governing 
function  for  the  variational  principle  is: 

N 


i-l 

A’- 1 


d(*) 


A' 


+  2^<u  ^  a  >a,+  2^<«/  .C 


(t)  ^UY 


il-l 

A’-l 


X'2 


(2.85) 


—  2  —2  ■¥  Boundary  terms 

Substituting  (70),  (83),  and  (84)  into  (85),  the  explicit  form  of  the  function  including 
the  boundary  and  the  discontinuity  conditions  is; 

N 


n = £{<  2r"'> 


(*)  o^k)Ak) 


Ak)  I  t^xAk^ ^ 


„(*) 


k-l 


+  <  -e*  +  £*  —  e*  5**^  k 


Cffiy^  yt 


otfiyb  yb 


+  <  +  +  2g*GjJ^  +  2Z^*^> 

+  <  +  g*  +  2g*  +  2x''^> 


,(*) 


+  <2€‘:J,-g*G'‘;33,3(2e;‘;;)+ro::'>  .. 


M) 


(/■)\ 


40^ 


(2.86) 
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+  <  Q1'  \  ^  -  g*  w  ;  >  +  g*  ( 2  <  ; )  >  , ,} 


UK 


JV-1 


j.  ^  WL<*)  <  j.  ^  1 

+  >j^*)->-<^  >g*<^33 


TV 


i-2 

AT-l 


+  2-^<«^a3  >  y«  +  <  0^33  >  ^)» 


i-1 


jpU> 


-2<v^  ,g*(r^3  >^, -2<w  ,^*0-33  >^, 

-2<V,  .-g*«r„3  >^,,,.-2<<f>^  >J^V)-2<>V  .-g*0-33 


+  £«  '’“^*^'^>“»  +  2*'rfV>s<.)+  -2<’V>S»> 

t-l  *  * 

+  <♦“^^^C',‘V“  +  2A'«r|P>^,+  <<^.g^C^V"-2♦"V>sy) 

+  <w“’,j<d‘‘’Q"’  +  2a.i)_)>^.,+  <Qf  .g^ci‘’w“-2«“\_)>^j) 

+£l<»".«H(c‘,*>y“)-  +  2iV;>V>s«)+<"“^*^<^*'’“*)'-2<’V>j<.. 

*-I 

+  <*"\*^(c','’<i)'  +  2j»''“v>4»>+<'<(-«*<<’«“’r-2<’n,)>^« 

+  <  w“’.«^(c‘5‘’e“’)’+22'.%)>j,„+  <e“’.s*((c'’’»'‘*’)’-2»'“’ii,)>^„) 
Here  we  have  used  the  notation  cr^*’  for  as  explicitly  stated  in  (52)  and  (55). 

Theorem:  The  Gateaux  differential  of  the  function  fl  defined  by  (86)  vanishes  if 

and  only  if  (70)  along  with  (83)  and  (84)  are  satisfied. 
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SisqI- 

Let  =  be  an  admissible  state 

corresponding  to  the  set  of  field  variables 

I  (*)  ^(*)\  _  J  .(*)  J*)  »/*)  Jt)  Jk)  ^U)  (k)  (k)  ^(k)  ik)  (i), 

{«  .<r  t  =  iv^  .c„3.!2„  .0-^3. 0-33} 

Gateaux  differential  of  the  governing  function  (86)  along  the  path  u,  provided  the 
limit  exists,  is 


8_ n(u)  =  .  g*  +  2  g*  ’  +  2 


n(k) 


(2.87) 


+  <^  ’8*Kps~^k\  >j,u. 

oP’S  ■'^ap  «  oPy6«  >8  6  "  oSyS'^yS  ;jU) 

+  <<*'.»•  <i,-<*Q‘.-*'‘’vr-/'*V“’  +  2j*c'*’  +  2z‘‘>>^„ 

+  < >y.> 

+  <  g*  +  2  f*  f'j* V  2  >  jjtt)  +  i  ~  y * V*’  > 

+  <  2<*3.f’Cl*’-^c'*’<.3y3<2<3)  >,(*)+  <  2<*i.g*Sl*’-g*G'*>„3^3(22‘;’) 
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N-l 


»(*) 


*-» 


a.  ^  J.  ^  ^  ^  ut  \ 

+  <  •J^0r3 


+  £f<  <'.-«*  <3*’ >^.,+  <  »“’.-|*<,  > 


*-2 


jjU) 


+  <  ,«} 


U)  «.-(*). 


A’-l 


+  2L'<^o3  >j,U)+  <0^33  > 


,(x)  ...Uo) 


,+a)  u)  *  u-ti) . 


o(t) 


+ <  <3’.w:’-K.*^“-K'"’>y..+ <  <’■«•»“’- W'’>^.,) 

-2  <  <.'’.8-<,'’ >,,-2  <  >^ 

-2  <  <"’.-*•  o-i"’  >y«-2  <  €''-‘;,g’<'T  >j^«-2  <  crjj"’ 

*-l  >  ' 

+  <  W‘".S^«"’t),-2#"t),)  >^„+  <  >^„ 

+  <  -e“’'l.  +  221‘’t).)>  J..)-  <  "’“'•S'5‘">1,>^ 

+  <tf.*’.S*<»'“’3).-2*"’31.)>^„+<e",sV‘>3,,>^,) 

+  £(  <i'l‘>.j^-(<’rT,,+2(;»“>rT,,)>  -  < v«>.g*(i5<“r„, > 

*-1  >' 

'2/  ^2/ 
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+  <  H‘“.*^(*'‘Vr,j-2(«<*V  V  >^„+  <  Af“;.ir(<'’)'T,,  >^., 


Similar  to  equations  (80)  through  (82),  the  following  relationships  hold 

<  >,.« = -  <  <  vr.s*??i>,>^., 

<  “’“'•s'si*’  <“’“*■**21;’’),  >j... 

<  ’1.  >j<« 

<  e“.i*w'“  >^)  =  -  <  q"  •«*«“’>/.)+  <  >^„ 

■where  =»  Boundary  of  U  internal  surfaces  in  Using  these  to  eliminate 

tL>  <.  3C,  «"1  Si  from  the  Gateaux  differential,  one  can  write 


*-i 


(2.88) 


^  -<*)  _*  D<*)  _*  n<*)  J*)  >. 

+  o&yl^yb-f^  ofiyh*^yh> 


O0y6  yb 
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+  <  2^3  .«•  e„  -  g*  G  (26^3)  > 

i(*)  ^/•-»-(*h  _*j.<*)  ...<*) 


+  <  S'‘’.**<2«"J)-g**“'-S*“'“’  >^ 

+2  zi<  <'.*•<>■:',“  >^,+ <  ?“’.>,  s*<'rf  >y.) + <  w'*’.«*4‘’ 

t-1 

+2  £{  <  \  ^  >  ,U)  +  <  -  g*  0-33  >  ,<*)  > 


jtf-i 


X-2 

A’-l 


,»  f  *  U)  .  ,  *  .«)  _*  U+l).  ,  ^:=+U)  _*.  U)  _*. 

>,U)+<®^33 


,(*) 


i(-l 


-2  <  >y2.)-2  <  >^Ar)-2  <  : 


A7 


*-l  ‘ 

**2 

+  <  w^*\  gH  -  t,^ + 0„  n„)> 

+  <  si;  \  g*( j  > 


4*) 


+2  £{  <  yi;\g^  -(<)'  +(0  »s<o 

*-i  ‘' 

+  <7/‘i,g^(v“’rT,,-(e>^)>^., 
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^4i 

+  <<!'‘',g^-(e“)'Ti,+(6jV>^.> 

^Si 

+  <2‘:’.«^(“'''’>'’l.-<*''’)'\)>s<;. 

Linearity  of  the  bilinear  mapping  <  ,  >  implies  that  vanishes  if  all  of 

equations  (70),  (83),  and  (84)  are  siilisfied.  Also,  nondegeneraleness  of  the  bilinear 
mapping  <  ,  >  implies  that  the  Gateaux  differential  in  equation  (88)  would  vanish, 
along  an  arbitrary  path  {iZ},  only  if  all  the  field  and  continuity  equations  and  the 
corresponding  boundary  conditions  are  satisfied. 


2.3.2  Extended  Variational  Formulation 

The  solution  of  (70)  must  belong  to  the  admissible  domain  of  the  operators  for  the 
functional  (86)  to  be  meaningful.  Using  (A.21)  along  with  (A.22)  relaxation  of  the 
order  of  differentiability  of  either  component  of  the  following  pairs  is  possible,  viz., 
(Q.  or  w),  or  ^„),  and  or  v^).  However,  relaxation  of  differentiability  could 
not  be  done  on  both  components  of  a  pair.  Some  possibilities  are  indicated  by  (80) 
throughd  (82).  Using  these  equations  to  eliminate  one  or  the  other  of  the  adjoint  pair 
of  operators  results  in  reduction  of  the  differentiability  requirement  on  some  of  the 
admissible  field  variables.  This  provides  a  basis  for  extension  of  the  variational  princi¬ 
ple  to  a  domain  where  the  differentiability  requirements  are  selectively  relaxed.  These 
extended  variational  formulations  also  provide  the  basis  for  certain  useful  specializations 
to  reduce  the  number  of  field  variables.  To  develop  these  extensions,  it  is  convenient 
to  rearrange  the  terms  in  (86)  to  write  the  governing  function  as: 


fl(u)  =  2  ’ ,  g*  *  +  t'"  > 


0  +  -g 


,(*) 


<- 1 
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+  <  w  ,g*7>3  +X  >  J 


U)  „U)  U)  D«)jU) 


+  aPy8V>j^*) 

+  <  > 

oifi  *  ^  (o>./5) 


+  <0^  ,g*M^p^-g*Q^  -k  V,  -y  <i>^^ 


X  ^  «)  *  »V'‘>  »  n<'^  »n<<> 


+  <  e* —  e*<t>'^^  > 


«)  * 


n 

.  ^  I  ^  X  ^ _*  \ 

■*■^^'^''0  ’  S*^o3  ^ 


-2<  vl‘\g*<r^3‘'  >,.-2  <  >,. 


^2<  >  ..,  +  2<  >  ...  +  2<  w'"'.g*<r-r  > 


(AT)  *  _+(Ar) 


+  r  <  vl^g*(  -  <>,  +  2  vl‘> - 


+  <^l''.f’(-'W‘“T]p  +  2y(!f'„’^7|p  <  Myp-np,g*(<f>V-^^V'^ 
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+  <  vv' ^  +  2  ^„ )  >  +  <  or  \  ,  g*(  w' ‘  -  2  .V '  > )  > 

St  •) 


U  ) ' 
6 


,(0. 


sU) 


(X) 


U  )’\'  ^  ) ' 


i(  =  l 


+  <w''>.8x-(e“’T)j+2(e“')-v>..,+  <e“\.s*((>»“’)'-2<*'‘ ’>')>,,.; 

•5c. 


->(*). 


^a) 


,(X) 


UK 


Ai  )\f 


STEP  1 


For  example,  eliminating  in  (89)  by  using  (80),  the  domain  of  N\‘il  is  extended 

from  C'  to  C,  and  the  function  has  the  form 

ft ,  (u)  =  2  £{  <  \  g»  4  \  g*  cr  4  z'^  S  ^ 

/t-1 

+  <  g*  4  } 

+D<vi*\-p^*4r’-/?‘*4r4 


*-1 


e  oP>6  >6  *  cpyb'^yb  ^  jfU) 


+  2  <  p*  —  p*  > 

jufO)  _»nO)  „u),u>  /UiZO 

<<f>,  -K  v^-J  0,  >^u, 


■*■  ■^  '^afl  ’  ^  ■^afl  <*^>6 ^  U*) 


,(*) 


+  <4‘\g‘Qr>F^*4<'4^,,, 


-h<Qr’.f*(2€r3)-g*<^r-g*vvr4^,,} 


(90) 
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A'  I 


+  L<<''a  ’«  *^03  >/*)+<  *^0  •^«*0-o3  >j,U>+<>V  .«*0-33 

<  I 


*-2 

A-I 


+  2-i<0-^3  •^'’a  -«*V^  >  +<0-33  ,g*W  -g*W  >  } 


Ai)  »  U+1) 


i-l 


-2<v^  .g*0-„3  >,,.-2  0^  .«*<^33 


+  2  <  vl;' \ g*  >  ^,,  +  2 <  crl'f '  > +  2 <  w‘-'  \ g* cr;^^ 

+  Z<  <vl^2g* <>,>  +<<>,. 2g^v;:>-o<:>)> 

/  -  I  «  '2 
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(A)  .  _-.(.V) 


+  <w*'’.8^-Q"7|^  +  2(3“’d.)>^„+  <eI,'’T)..S^t>’“>-2»'‘>)>  J 


+  Z<  <  vr  >^.,H-  <  <’.,.2*^ /,*'-<>)  >^„ 


k-i 

lU)  .*/  MM 


+  <-'‘’.^-e-l‘’’>,+2C’i.)>^..+  <el‘S..s^w''‘'-2.»'“>)>,.,i 

■^5/ 
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STEP  2 

Eliminating  in  (90)  using  (81)  leads  to 

N 


*-i 


,(/.) 


+  < 


N 


w'*‘.s*J^"  +  x'*’>  ,.,1 


+  £t<v“>,-f“’vl‘’-/i''V" 


Jk) 


t-\ 


+  <e 


afi’  *  ofiyhS  ^yl.  6  " 


+  <  >,u 


nU) 


+  <2^3^f*Ql*^-^G<^3y3(2<3^)>^„, 

+  <  q1*^  .  g*(2  <*J)  -  g*  <  -  g*  } 

\-\ 

+  L^<v„  .g*o-„3  >j^a)+<<^o  ’'<^‘^03  >/*)+< 

*-l 

+  £{<  vL*\-g*<r:r>^,+  <  >v"\-g*a- 


(91) 


w 


u)  *  -(/).  1 

-^°'33>„U)J 


*-2 

N-l 


*-l 

^  <l)  *  -(*)  ^  1  ^  <l)  *  ^ 

-2<v„  .g*a„3  >  -2<w  ,g*a,,  > 


30 


.  (A)  *  -'(A).  .  -%^jl(A)  .  *_‘(A)^  .I-.,'-'') 

+  2<v^  ,^0-^3  >^,^,,  +  2<0^  -fvg  <^o3  >j^ia.  +  2<w  .  +  g*o-33 


+  >^„+  <  >5(*) 


a)  ^UK 


+  <  ».''>,j^-c",,^  +  2a"v>j.>+ <e“\.S^“’'*’-2«“’)>j.,  ) 


+  £(  <  vl'',2s-r>,  >^„,+  <  a;.>,.2s-(  V“’-<'’)>^,, 


+  <<",2s*J»'“;ti,  >^„,+  <  m1>,,2s^ <*’-<'’)  >^., 


STEP  3 

Eliminating  j2^*^  in  (91)  using  (82)  results  in: 

fi3(u)  =  2 £{<  >,(«+  <  > 


o  Q  ^  a 


oPy6^y6  opy6^y6  ^ /t) 


+  2  <  g*  —  £*  > 


’  Of  ^O  j^ik) 


■*■'^'^0^’  00^6  Sfi  a^ye^e  ^  nft) 


(92) 


+ 2  <  .  g*  ( 2  )  -  g*  w|>  ^,, , } 


N-l 


t-l 

N 


..(*)  _»  _-<*)  ^  _L  ^  ...(<) _ *  ^  I 

33  r<*> 


.  T'/^  UJ  _*-«)-  ^  »_ 


k-2 

N-\ 


LI  Hi)  ^  U)  ,  i  "•l>,  ^  ,  ,  ‘Ui  *  UJ  i 

<<0'«3  +^^*^0  ~S*K.  +  >,.U.+  <°'33  ’g  ~S*^'  >^,i 


,(/t)  *  (/•!) 


(.')  *  U)  *  (ii'l) 


/-I 


-2<  >^,  -2  <  w >^, 

+ £{  <  >^*)+ <  0,.2g*(  v;;^-ol*') 

*-1  *  * 

+  <*“’,2** *2^,  >^^.,+  <  Af">'.),.2*^*|;'-0“’)  >^„ 

+  <w"’.2j*fi“D.>^.,+ <e“N..2s^»’“'-«''’)>^,) 


+  £{  <  vr.2g*i^'l>,  >  .,+  <  Nl>,.2g*(  V--0-)'  > 


«). 


*-i 


V*)  '  ^  '’aP  ’'^’*'6  '  "o  'o  ^ 

•’a 


+  <«"^2y•A■“’7,J  >^„+ <<;»!, .2r(«l‘’-<’)' 


+  <w*.2s"fi-:;’D.>^,+  <erT,..28Hw"'-*"')-  >^) 


2.3.3  Some  Specialization 

The  function  defined  by  (92)  has  no  differentiabiJity  constraints  on  the  stress 
resultants  N This  function  also  leads  to  certain  useful  specializations  by 
identically  satisfying  some  of  the  field  equations.  Assuming  that  in  (92),  continuity  of 
tv'*’  is  identically  satisfied  i.e.  normal  displacement  is  restricted  to  be  constant  through 
the  thickness,  i.en  w'*’  =  w  for  all  k  (67), 


O  o 


+  afiyS^yi  afiyb^^yb  ^ /H) 

+  2  <  e*  —  e*  > 


+  ,-i?  -/ 


+  aPyb^yb  ^  afiyb'^yb  ^  gik) 

+  <2<^],-g»G^'>^3^3(26<;>)>^,, 

+  2  <  e“'  ,*•  C2  «“>  -  j*  «"-  *•  w“>  >  ^„ ) 


<*)  . 
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+  £«  >^.,+  <  .r,;  >^.,) 


+  L^<^.3  'S  ^’o  -g*\  +  >.,u.+ -8*^  >,.2)5 


,(2)  *  (/*!) 


3.^ 


-2<v^  ,g*<r^3  >^,-2<w  ,g*(T^^  >^, 


■*■2'*'  '’a  ’^^aS  ^  „<A’)  ■*■  2  N  ^  ^  ai  ^  j,(A’)  ^ ’  “*"  S  ^33  ^  i>' '  ' 


+  ZI  >^0+  <A'"»’)j- 


it-1 


4^' 


.XA) 

^4 


+  <0r.2g» 

^3  ‘^4 

+  <  'v"\2g*0r%„  >  ..+  <Q:;S..2g*(>v<‘^-v^‘^')  >  } 

‘^S  '^6 

+  h  <  <  ’  ’  2  g*  T,^  >  ^.3  +  <  n;:;  .  2g*(  >  -  0;: ' )'  > 

•>».  'n. 


/-I 


+ «^::\2^o,  <  MST),.2g*(0‘:^-<v  > 


^  w^.--6  '^y'o,  -  ^k) 

•  4/ 


If  the  kinematic  relations  (59)  through  (61)  are  identically  satisfied,  equation  (93)  leads 
to: 

N 


*-I 


,(*) 


N 


+  J-(<v“>,-f‘"v<‘>-/i'*V‘>> 

Of  o  ~a 


,(*) 


ft- 1 


olfiyh  yb 


+  <*".-*“'v">-/'V'’>  ,„ 

^or  '  a  ^or  p\K) 


(/)  f-.j'i) 


+  <  2<3\-g*G-„3^3(2e';')>^,,,+  < 


(94) 
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+ 1<< 


+ D<  .-g*<r  >,(*>+  <  >,u,} 


+  2  <  .  g*  <T:<f  >  +  2  <  <f>^:\t,  s* crl'f  >  ,.,V.  +  2  <  >v'^’^ .  +  g*  cr;'f  >  ^,,, 


+  £{  <  v;  \  2  g*  r,^  >  ,. ,  +  <  1, .  2^  ^ 


+  <  w<*\2g*0l*^,  >^«+  <Gl'^««2g^>v^‘^-t^‘*^)>^„  } 


+  D  <  >^*)+  <  ^^^^.2g*(  v«>-i>^>)' 

*-i  *'  ** 

^3l  * 

+  <  w<‘\2g*0'l%„  >^o+  <Gr^„.2g*(w^'^-^<‘>)'  } 

If  the  kinematic  boundary  conditions  in  (83)  and  (84)  arc  identically  satisfied,  equation 
(94)  reduces  to 

n,(u)  =  2£{<  v<;\g*y*4y<*’>^,,+ 

*-i 

+  <>v‘*\g*F‘;^x^*^>  ,} 


I  T'  I  ^  *'* '  o« )  «  )  tjU )  j « )  V,. 

+  >(<v  ,  —  P  V  —R  9„>( 

MM#  O  <»  ~  ^  iP' 
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^  o/fiyb  yb  ^  afiyb  yb  nU  ) 


+  <  ,  -  i?*' '  ’  0^'’  >  , . , 

^a  ’  a  ^or  j^U) 


+  <  >,U> 


+  <  2€S.-g*G^'U3(2€‘;30  >  +  <w.-P''"w>  } 


+  L^<''o  *^^<^03  >/*)+<<^o  •'/«  O’o.l  >,,<*.+  <"’  •«  ^33  >,,w 


U)  ,  ..  +(/ )  . 


(O 


+  L^<''o  *-^^03  >„U)+<'^  •-?^0'33>„u.> 


+  L^<0^«3  •«*'’«.  +^«*^a  >  „U)+<  0^33 


,(*)  »  (/HI) 


-2<  vl‘\r  o-i*^  >^-2  <  >^, 

+  2<  v!f \g* >  .A,)  +  2<  C^h ^ >  ,CA.)  +  2 <  »v‘"’\  +  >  ,,, 


+  2£{  <  >  ,*,+  >  ,*)+  <  } 


+  2  £{  <  >5(*)+  <<'.«* O,  >,<3.+  <  '^“^2g*<5-l'S„  } 


^3-°^33  "®^33 

Then,  if  the  constitutive  relations  (62)  through  (64)  are  identically  satisfied,  Sl^u) 
specializes  to 

0,(1.)  =  2£«vf  .**y'’  +  I^'>>  ,+  <*f  .g*G“>  +  z'"> 


+  <w.«'f"’  +  x“’>  .„! 
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N 


+  £'  -  <  <  •  S*  C  <1  ■  s* 

/'J 


+  <9o,  .-f*(2„  v^  -/  >^a)+<w,-P  w> 


/*) 


(96) 


+  L^<  '^1*  Vo-^a'’  >,<«+  <  >-(*)+  <  „u)J 


/‘I 

N 


+  L<  <  ’'o  ’ '  -  3  '  >  „U.  +  <  '^  •  -  g*  0-33  >  p(*)  } 


1-2 

A  1 


+  2-^<°^oJ  ’^''o  +‘*^*<^0  -g^'’«  >^U)} 

i-1 


(/  )  *  .«  >  1 ) 


-2<v^  .g*<r^3  >,.  +  2<v^  .g*o-^3 
+  2  <  *1"’ . S*  <,"’ >  ^„,  +  2  <  w .s*  J>3  > , 


+ <  ^'1*  V  lip  >  ^*)  +  <  g*  ■n^  >  s(*)  +  <>v .  g*  0,  -o,  >  ^*) } 


i-l 

N 


4*> 


4*>^ 


+  2£(  <  .  g*  Kl  >  s(*)  +  «f>f  •  g*  ^ t2  Tip  >  X*)  +  <  'V .  g*  ii„  >  ,*, } 


X-I 


If  the  interlaminar  displacement  continuity,  (66),  is  identically  satisfied,  (96)  reduces  to: 


N 


n.w  =  2£«v«.**i^«+i'“>>  +  <«;;',j.c«+z“’> 


*-l 


ffk) 


+  <w.(r*J='‘’+2f">^„) 


AT 


+  £{-  <  w'«.j.v"  ,>^.,-<  ,> 


,(*) 


t~l 


-  <  G» ’.g*''’.o>y*.+  < 

+  <0  ,  —  gQ  —K  V  — /  ©  >,, 

’  O  *^o  O  ^o  y*) 
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(97) 


-<w,P 

(1)  *  -<1).  .  (-W)  *  '(.V) . 

-2<v^  .«*o-„3  >yA) 

+  2<«“'.^S«<r;,''’>  ,„+2<  >, 


+2Z(  <»“.*•  >j,.,+ >,,.,+ <  «-g‘6,x  >.. 

t-1  •  3  S 


+  2^<  <  >  .*)+  < 

/«!  I'  3'  S' 

Here,  v^*'  is  no  longer  independent  for  k=  1,  2 .  n.  I'urthermore,  if  siirlace  sliear 

tractions,  surface  shear  couples,  body  forces,  and  body  couples  are  neglected,  i.e. 

=  =  =  =  =  0  then  (97)  reduces  to; 

=  2 £{< >  ,„+  <  >  } 


+  ^{  <  ,  g*  ,  g*  yO 


•  O  >C.q,  O  j^ik) 


-<w,F*’w>  ,,,} 


+  2<w,g*P.>, 


+2£{  <  >^*,+ >^o+  <  >^,,} 

*-i  ‘  *  * 


+  221{  <  >5t«+  <<"•?* >5(*)+  <  >  ,.*)> 


38 


Explicitly,  replacing,  v*'  using  the  continuity  of  the  displacements  across  the  interface, 
by  v”’  and  fl,  is  written  as: 


(/)  vrU) 


At)  r,<t)  , 


£«  Z',c>  > 


(X)  ^  ,  ,(/) 


.  -  .  ( 1 )  r^i)  JO  ^  ^  Jt>  JO^  ,  A')  ^ 

+  2<v  ,R  <b  >  ,,.  +  2<0  ,K  >  t  d>  >  ,, 

n  ~o  nit:)  l~a  -At 


t  1 

U  )  ,Xi )  ^  .  ,(/) 


o  o  n't)  a  ’  1  a  At 


t-l 

(1)  o(X)  ^  ^  ,(/) 


I  ^  ^  )  I  )  -iw  \ 

+  <9„  .rQ»  +-'  >  j*)  / 


+  2<w.g*F,>, 


+2£{<  >, 


+  <  >At>+  <  } 


If  inertia  terms  are  ignored  i<.  for  the  static  problem,  it  is  not  necessary  to  take 
convolution.  Then,  the  governing  function  is  the  functional: 


t  I 

,«)  .  A') 


1  /“)  =  -  ’  S*  •  «*  L  Ko  >  r.' 


(1(K)) 


+  <  ^a^i’^'^(>.W>j,u.+  <C>..  ^8  „<n+  <^..  ’8  Q,.  >^a,  } 

+  2<w, g*P^>^ 

+  2£«  r A/l>,  > , 


<-i 


•  '^‘/-^a’S  •’aP'lfi-'^t) 

*  i»l  ’ 


+  <  >^,,+  <  'V.f*0„T)„>^,U,  } 


/  I 


+ 22:{  <  V- ,  +  <  £ ,  r:;  ^ 

X  I  '  /  I 

+  <<t>\'\g*  /fir'l'p  ,  +  <  H- .  g*  0',.  T1„  >  ,  I 


.{<^ ) 

1; 


If  the  physical  problem  does  not  have  line  loads  or  couples  applied  to  surfnces  in 
interior,  and  ^  vanish.  This  gives  the  governing  functional  as 

N  X-l 


*-1  (-1 


+  2<w,g*Pj>j^  (101) 

This  specialization,  dropping  the  convolution  with  g(t),  was  used  to  set  up  the  finite 
element  approximation  discussed  in  the  next  section. 
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Section  III 


FINITE  ELEMENT  FORMULATION 


3.1  INTRODUCTION 

I'he  I'inite  element  method  subdivides  a  given  region  R  in  an  n-dimensional 
euclidean  space  into  a  number  ol  disjoint  open  subregion  (elements)  R' ,e  =  \ 
such  that 

R=  Vim  M  W  (102) 

— ijii 

Here  a  superscripted  bar  over  a  quantity  denotes  its  'closure'  i.e. 

S"  =  U  9^"  (103) 

where  QR'  is  the  boundary  of  R'.  Disjointness  of  these  elements  implies 

R‘'P|r'  =  0  if  €9tf  (104) 

A  set  of  nodal  points  in  O  defines  the  geometry  of  the  elements.  In  the  following 
section  a  discretization  of  the  domain  by  the  finite  element  method  is  presented.  The 
formulation  of  the  finite  element  is  based  on  the  variational  principle  governed  by  the 
specialized  functional  fl,j. 

3.2  FINITE  ELEMENT  DISCRETIZATION 

Let  the  field  variables  at  any  point  within  an  element  be  represented  by: 

op  V  p  Q 

</>‘^'(xJ  =  {//^(a-  )}(«!>"’}  (105) 

op  ^  p  o 

MxJ  =  i//'(x..)}{W} 

P  '•  P 
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Here  {/7,} ,  are  sets  of  interpolation  functions  relating  tlie  \alues  of  the 

corresponding  variable  at  the  nodal  points  to  an  arbitrary  jxiint  within  the  element  and 
«*’},  and  are  vectors  of  the  apriori  unknown  values  of  ,  w  at 

the  nodal  points.  Substitution  of  (105)  into  the  governing  functional  (101)  yields  the 
spatially  discretized  functional  as: 

m 

n  =  L  n, 

(  I 

where,  for  N  layers  in  llie  laminate: 

A’ 

'  1  p- 

I  I 

tt 

if 

+  2  r {V^'YiDH  ][B^^){DHf{^^*)dR" 

Jo  V  <t>  y 

if 

t-i 

+  2/*  (106) 

If 

If 

{DHJ  (W) 

+  f  <{wf  A*^'Y  >  [C^'){{DHf,{HY]  dR'  } 

i  “  {//,} 

+  2  J{wf  {HjP^dR' 

if 


Here  the  symbol  D  associated  with  H  denotes  appropriate  differentiations;  [/A**’],  [^‘’] 
and  [D“T  represent  constitutive  relations;  and 


=  r  {DH  ]lA^‘){DJl  dR‘ 

=  f  {DH )[/^^{Dn dR" 
pT 

=  /  {DH^}[A^%DH/ dR^ 
S' 

=  f  {DH  )[B^%DH  f  dR' 

Vtp/t  I  V 

s/‘ 

J  {D/JjlB")iD//j'dR' 


s 


f  {DHjUf)iDH/dR‘- 

=  r  {DHj[C^‘^){DHfdRr 
jf*) 

R  =  f  {//  }P,dR‘ 

I.  /'  J  '-3 


where 


[Df//  = 


ID///  = 


Ih^ 

0  1 

.■A 

0 

//^ 

V.> 

//^ 

‘'.J' 

1%1 

Ih^ 

0 

<p.> 

0 

K 

<p-> 

//^ 

0,V 
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{dhJ  = 


h],  0 
♦ 

lo  h] 


{H)  =  {//.,} 


Using  the  above  definitions,  equation  (106)  can  be  written  as; 

.V  I  I 

c  a  ttci'*  y  t\  *■  \ip,r  t  > 

/  I  /I 


/-  I  i- 1 

(-1 

+  +  {Wf  {W} 

+  (w)  +  {wf 

+  }  (107) 

+  2{wfR^,^ 

Gateaux  differential  of  the  governing  functional  (107),  denoting  by  {FJ,  {$^},  {VV|  the 
path  of  variations  in  {V^} ,  {<I>^} ,  {w}  respectively,  is: 


1-1 


I- 1 


i~\ 


o  >  ipn  )  <1  y 
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(108) 


<  1 


/  1 


I  1 


/  1 


+ {55;jy  [<^f  w +{m''[<;^]{<i>5;^ 


+{^<;y[OK^ } 


+  2  R 


"7' 


The  vanishing  of  the  Gateaux  differential  S^ft  (108)  for  arbitrary  {FJ ,  {<I>J ,  {W} 
yields  the  following  equations. 


C-  1 


/  1 


AO 


7-1 


1-1 


=  0 


(109) 


m  N 


£1  £(-1*'.^ 


»-l  7-1 


<-l 


(no) 


N 


7-1 


7-1 


XMCfiv'JV 


7-1 


7-1 


7-1 


1-1 


£<,iC)£M'»'^-£‘<iO(<’> »  = « 


I- I 


and 


£ 

r-  1 


£l-U 


(7) 


AO, 


7-1 


=  0 


(111) 


Where  the  ^  denotes  the  ’’direct  stiffness”  addition  of  contributions  of  all  elements. 

r*l 

Equations  (109)  through  (ill)  can  be  written  collectively  in  the  following  matrix 
from. 


[K]{U)  =  {R] 


Here 


i/f]  = 

c-  1 
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{/e}  = 


where  n  is  the  total  number  of  nodal  points  in  the  system  and  we  identif y: 


^.2 

^.3 

/f,,  ..  .. 

••  -  ^,.V 

^..VU 

^  l,A  *2 

^22 

^23 

^24  ••  •• 

^2^*1 

••  ••  ^2.A- 

^2.A.. 

^2,A>2 

^33 

^34  "  “ 

••  ••  ^3.;v 

^i.N*2 

^44 

••  ”  -^4JV 

y 

If 

■^4,A?+2 

symmetric 


K  K  K 

^N.N  "A^,A’t2 

If  If 

A’+1,A*1  A-H,/Vt2 


where 


=  Zi<!; 


jr„  =  o 


1-2 


46 


,zio 

/  / 

iio 

l~N-l 

*-o-., =iC’'+v-.ito 


t^N 


*,,V.3  =  lO 


^^33  =  LK'J 


/  1 


*33=i<';i 


*■34  =  lo 

^33..  = 

^3,. = 

^3.».,  = 

^3_V.3  =  1<«1 


N 


^33  =  io+iO+':ZfC.J 

hi 

^3-=',lO  +  '.‘3ZlO 

/-3 

/-/ 

^.»=',(C’'  +  '.'.V3  Z'O 


/-.V  I 


^3,.V..  +  ,£[<;,] 
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^3JV.2 

*42..  = 

l-i 

*4.A'  =  '2l*i^'’'  +  '2',V  2  Z  lO 

/- Y  I 

*4A.,  =  ‘2l*!i,.'’l  +  '..‘v  ,ZO 

/  A 

*4J.-.2  =  ‘.lO 

/-» 

*„,j,  z  lO 

/•AT- 1 

*AJ.  +  Z 

/-JV-1 

*AJ4..  ='A-2l*Si'*l  +  '»-2J»-.ZlO 

*AJ4+2  *  *W-2^**4J 


.V 


/-JV 


/  .V 


K 


N^  i,A'+a 


= '.V 
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and 


n;)^  =  <v;.w.*;,0= . 


U),= 


'-p\ 


K) 


lo 


where  the  subscript  j  denotes  the  jth  node. 


(112) 


3.3  SOLUTION  PROCESS 

Since  the  stiffness  matrix  is  symmetric  and  banded,  only  the  above-diagonal  terms 
in  the  global  stiffness  matrix  are  stored  in  the  form  of  a  rectangular  matrix.  The 
first  column  of  this  matrix  is  the  diagonal  of  the  original  matrix.  The  load  vector  is 
set  up  according  to  (112).  After  the  global  stiffness  matrix  and  load  vectors  are 
as,sembled,  the  system  of  simultaneous  equations  is  solved  for  the  displacement  degrees 
of  freedom  using  the  Gaussian  elimination  process.  The  solution  for  nodal 
displacements  (S')  of  an  element  along  with  [5]  matrix  which  is  the  element 
St  ram -displacement  relation  are  used  to  calculate  the  stresses  at  the  desired  points  of 
the  element.  i.e., 

n 

e  =  BB"  BB 

I  I 

/-! 

where  n  represents  number  of  nodal  points  in  the  element, 
a  =  QBB' 
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and  [Q]  contains  the  elastic  constants  of  the  specific  lamina.  In  (106)  the  specific  [Ii\ 
is  denoted  by  [DllJ 
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Section  IV 


VERIFICATION  AND  APPUCATIONS 

4.1  PROBLEM  DESCRIPTION 

The  theory  discussed  in  Sections  II  and  III  was  implemented  in  a  computer 
program.  The  computer  program  was  written  in  I'ortran  and  designed  to  run  on  the 
IBM3()81  main-rranie. 

r.ight-n(xled  .serendipity,  nine-noded  I.agrangian,  and  four-noded  Lagrangian  quadri¬ 
lateral  isoparametric  plate  elements  were  used  (Fig.  2).  The  nodal  degrees  of  freedom 
(the  number  of  unknowns  associated  with  each  nodal  point)  depend  upon  the  number 
of  layers  in  the  laminate.  If  the  number  of  layers  is  m,  then  the  degree  of  freedom 
is  (2m+3)  for  each  nodal  point.  This  consists  of  three  displacement  components  (u,v,w) 
for  a  reference  point  on  the  'through  the  thickness'  nodal  point  and  rotations  (0,,  0^) 
of  each  layer  to  describe  completely  the  deformed  geometry.  Two  versions  of  the  pro¬ 
gram  were  prepared.  The  first  was  an  ’’incore”  program  wherein  all  the  matrices 
were  generated  and  stored  in  core  and  the  solution  process  did  not  require  auxiliary 
storage.  This  version  was  adequate  for  solution  of  example  problems  and  verification 
of  the  general  approach.  Later,  a  version  was  developed  where  the  algebraic  equations 
were  assembled  and  stored  in  blocks.  This  version  was  used  for  solution  of  a 
22-laminae  problem  using  the  CJRAY-XMP  computer  at  the  Ohio  Supercomputer  Center. 
This  implementation  was  verified  by  application  to  several  example  problems.  These 
examples  included: 

1.  A  simply  supported  square  sandwich  plate  made  of  isotropic/orthotropic  surface 

layers.  The  core  was  assumed  to  have  a  finite  shear  modulus  but  zero  elastic 
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1.  Element  04:  (a)Rectangular  parent  (b) Isoparametric  counterpart 


Figure  2:  Four,  eight,  and  nine-noded  quadrilateral  element  used. 
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modulus  in  extension. 


2.  Free-edge  delamination  specimens 

a.  Angle-ply  [±45], 

b.  Cross-ply  [0/90]^ 

The  results  were  compared  with  available  solutions.  The  effect  upon  the  accuracy  of 
the  results  of  mesh  refinement  and  sublayers  subdivision  was  studied.  The  code  was 
used  to  determine  stresses  and  del'ormations  in  a  multi-ply  free-edge  delamination  speci¬ 
men. 


4.1.1  Analysis  of  Plates 

The  multilayer  plate  I'inite  element  methcxJ  was  used  to  analyze  a  three  layer, 
square,  simply  supported  sandwich  plate  uniformly  loaded  in  the  transverse  direction. 
The  geometrical  and  material  properties  of  the  plates  were  (the  top  and  bottom  layers 
are  denoted  by  subscripts  1  and  3,  the  core  by  2).  This  example  was  the  same  as 
used  by  (5] 

Plate  dimensions; 

Length  of  each  side  =  10  inches. 

Thickness  tj  and  t^of  surface  layers  =  0.028  in 

Thickness  of  core  t,  =  0.75  in 

Material  properties 

a.  Stiff  layers  are  isotropic  elastic. 

£,=£3  =  10’  Ib/in^ 

V,  =  0.3 

=  3x10''  Ib/in^ 

b.  Stiff  layers  are  orthotropic  elastic. 
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\ 


,  =  /;’  =  Ax\()*'lb/iiC 

y  •  >  ’ 

G  ,  =  G  ,  1.875j:1()‘’  lb/in~ 

v>’l  lyJ 

1/  ,  =  V  =  0.3 

A 1 

G^.2  =  3x10''  W/in^ 

G  ,  =  1.2x10"  Ib/in^ 

yr.2 

lading 

^=1  Ib/in^ 

Due  to  double  symmetry,  only  one  quarter  of  the  plate  was  analyzed  in  each 
(I'ig.  3).  When  eight-  and  nine-noded  elements  were  used  four  mesh  discreiiza 
viz.,  (1x1),  (2X2),  (4.\4),  and  (8X8)  were  used.  In  the  case  of  four-noded  element, 
discretization  for  the  quarter  of  the  plate  was  extended  to  a  finer  (16x16)  mesh. 


case 

■ons 

the 


.S4 
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4.1.2  Analysis  of  Four-Ply  Free-Edge  Delamination  Specimens 

The  free-edge  delamination  specimen  was  treated  as  a  special  case  of  a  plate 
rig.(4).  The  boundary  conditions  considered  were:  fixed-fixed  (longitudinal  displacement 
specified)  at  the  longitudinal  ends  and  free-free  at  the  transverse  edges.  The  results 
for  four-ply  symmetric  laminate  with  stacking  [±45],  and  [0/90]^  were  compared  with 
Pagano's  solution  [10]  which  is  based  on  a  generalization  of  Reissner's  theory.  The 
comparison  covered  only  cr,  and  because  Pagano's  solutions  were  available  only 

for  these  quantities. 

The  dimensions  were  specified  as  a'b^  1  and  b  =  8h,  where  a,  b,  and  h  are 
length,  width,  and  thickness  o(  each  lamina,  resjrectivelv.  For  proper  comparison,  the 
material  properties  assumed  were  the  same  as  in  [lO]. 

=  20x10**  psi 
E22  =  £33  =  2.1x10*  psi 
G^2  *  G,3  »  G23  =  .85x10*  psi 
•'i2  =  *'|3  =  ‘'23=«-21 

Following  Pagano,  N«2  indicates  each  lamina  is  treated  as  a  single  layer,  and  N=6 
indicates  that  each  lamina  of  thickness  h  is  modeled  by  three  sublayers  with 
thicknesses  of  h/3.  For  N=10  each  lamina  with  thickness  h  is  subdivided  into  five 
sublayers.  In  going  from  three  sublayers  to  five,  only  one  of  the  sublayers  was 
subdivided  into  three  new  sublayers  with  thickness  of  h/9. 

Two  different  types  of  mesh  refinement  were  considered.  In  both  types, 
refinement  was  carried  out  in  a  manner  such  that  the  nodal  points  of  the  previous, 
coarser  mesh  were  a  subset  of  the  finer  mesh.  Refinement  associated  with  longitudinal 
direction  was  performed  by  dividing  the  domain  in  the  longitudinal  direction  into  the 
stated  number  of  equal  length  elements.  In  the  case  of  the  transverse  direction,  the 
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refinement  was  in  the  form  shown  in  Figs.(5)  and  (6),  to  minimize  the  required 
storage.  In  one  sequence  of  refinements,  shown  in  Fig.(5),  the  width  of  the  specimen 
was  first  divided  into  five  equal-width  elements.  Only  the  edge  elements  were  refined 
further,  each  time  the  edge  elements  being  subdivided  into  three  elements.  This  process 
was  continued.  In  the  other  sequence  illustrated  by  Fig.(6),  an  edge  element  was  u.sed 
such  that  the  (y/b)  ratio  for  its  center  was  0.995.  The  remaining  interior  domain  was 
discretized  into  three  strips.  Refinement  over  the  thickness  was  done  on  the  central 
sublayer  while  changing  from  N=6  to  N=10,  Fig.(7c)  in  order  to  get  good  estimates  for 
T,,  and  tr,,.  On  the  other  hand,  to  evaluate  r^.  at  the  interface,  discretization  was 
performed  on  the  sublayer  next  to  the  interface,  Fig.(7d).  Because  of  symmetric  stacking 
sequence,  only  two  layers  were  used  in  the  analysis. 

4.1.3  Results  of  the  Analysis 
a.  Rectangular  Plates 

Figs,(8)  through  (11)  show  a  comparison  of  the  approximate  solution  with  that  of 
the  sandwich  plate  theory  as  presented  in  [5].  Fig.(8)  shows  the  convergence  of  the 

central  deflection  of  a  sandwich  plate  with  isotropic  layers.  Results  using  the  Q8  and 
the  Q9  elements  converge  to  the  exact  solution  more  rapidly  than  those  from  the  Q4 
element.  The  results  obtained  for  one  element,  in  the  case  of  09  element,  or  four 
elements  in  the  case  of  Q8  elements,  are  superior  to  results  obtained  from  256  (16x16) 
Q4  elements.  Fig .(9)  shows  that,  for  a  given  error  limit,  the  amount  of  time  needed 
for  08  and  09  elements  is  considerably  smaller  than  for  04  elements.  Results  for  a 
sandwich  plate  with  orthotropic  layers  are  presented  in  Figs.(lO)  and  (11).  The  same 
observation,  can  be  made  as  in  the  case  of  isotropic  plates.  The  numerical  results  for 
isotropic  and  orthotropic  cases  are  presented  in  I'ables  (l)  and  (2),  respectively. 

The  program  was  also  u.sed  to  solve  the  isotropic  problem  where  the  shear 
correction  factor  was  ignored.  The  maximum  central  deflection  for  an  8x8  mesh,  using 
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(a)  Five  subdivision  across  the  width 


(b)  Nine  subdivision  across  the  width 


(c)  Thirteen  subdivision  across  the  width 


(etc.) 


Figure  5:  Sequence  of  mesh  refinement  across  the  width. 


60 


61 


MRXIMUM  DEFLECTION  xlO'^ 

..00  0.02  0.04  0.06  0.08 


X  Q9  ELEMENT 
+  08  ELEMENT 
A  04  ELEMENT 
O  SERIES  SOLUTION 


Figure  8:  Number  of  subdomains  along  each  side  vs  central  deflection  of  a 

square  isotropic  simply  supported  plate  with  No.  of  elements  >= 

N\ 


-f  Q9  ELEMENT 
A  08  ELEMENT 


O  04  ELEMENT 


Figure  9;  CPL  time  vs  error  in  central  de 
isotropic  layers. 
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I - ^ - Y 

0.00  40.00  50.00 

^P/23  (CRRT) 


with  mesh  refinement- 


MAXIMUM  DEFLECTION  kIO'" 

.00  0.04  0.08  0.12  0.16 


X  Q9  ELEMENT 
+  Q8  ELEMENT 
A  Q4  ELEMENT 
O  SERIES  SOLUTION 


Figure  JO.  Central  deflection  of  a  square  orthotropic  simply  supported 
sandwich  plate. 
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PERCENT  ERROR  OF  DEFLECTION 
.00  20.00  HO-00  60.00  80.00  100.00 


+  Q9  ELEMENT 
A  Q8  ELEMENT 
O  04  ELEMENT 


Figure  1  h  CPU  time  vs  error  in  central  deflection  with  mesh  refinement- 
orthotropic  layers. 


65 


Tabu  1:  Numerical  values  for  maximum  deflection  and  the  corresponding 

CPU  time  for  different  types  of  elementsXisotropic  case) 


um  deflection 


0.000064402 


0.00036986 


0.000S968S 


0.000699S9 


0.00072964 


0.0006189 


0.00073851 


0.00073997 


0.00074005 


0.00073539 


0.00074087 


0.00074008 


0.00074006 


0.00074 


Number  of 
Ekmenis 

Element  Type 

CPU  (Sec) 

1x1 

Q4 

0.046 

2x2 

Q4 

0.154 

4x4 

Q4 

0.605 

8x8 

Q4 

2.918 

16x16 

Q4 

19.897 

1x1 

Q8 

0.232 

2x2 

08 

0.934 

4x4 

08 

4.492 

8x8 

08 

29.377 

1x1 

09 

0.287 

2x2 

09 

1.221 

4x4 

09 

6.516 

8x8 

09 

50.446 

Series  ^ 

Solution 
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Table  2:  Numerical  values  for  maximum  deflection  and  the  corresponding 

CPU  time  for  different  types  of  elements.(C>rthotropic  case) 


Number  of 
Elements 

Element  Type 

CPU  (Sec) 

Maximum  deflection 

1x1 

Q4 

0.046 

0.00026536 

2x2 

Q4 

0.146 

0.00066988 

4x4 

Q4 

0.591 

0.0010218 

8x8 

Q4 

2.851 

0.0011671 

16x16 

Q4 

19.071 

0.0012078 

1x1 

Q8 

0.228 

0.0010439 

2x2 

Q8 

0.954 

0.0012194 

4x4 

Q8 

4.497 

0.0012215 

8x8 

Q8 

29.467 

0.0012216 

1x1 

Q9 

0.286 

0.0012140 

2x2 

Q9 

1.224 

0.0012224 

4x4 

Q9 

6.545 

0.0012216 

8x8 

09 

50.267 

0.0012216 

Series  Solution 


0.00123 


08  element,  was  =  ().(K)()43601.  I'lie  series  solution  for  same  problem  gave  almost 

identical  value.  =  0,(XX)436()5.  Therefore,  maximum  deflection  due  to  shear  for 

this  plate  would  be: 

W-  =  0.000740  -  0.0004360  =  0.000304 

^max 

b.  Four-Ply  Free-Edge  Delamination  Specimen 

Values  of  oTj,  presented  in  the  subsequent  sections  were  determined  based  on  the 
constitutive  equations  (.36)  and  (37).  Determination  of  these  stress  components  is  also 
possible  through  the  use  of  the  equilibrium  equations.  However,  inplane  strain  is 
linear  over  an  element  and  discontinuous  at  the  ntxlal  points,  resulting  in  a 
discontinuous  set  of  inplane  stresses.  Furthermore,  the  numerical  evaluation  of  O’  , 
needs  one  numerical  differentiation,  and  tr^j  needs  two  numerical  differenti;  ion 
Therefore,  the  use  of  numerical  differentiation  would  not  result  in  a  satisfactory 
estimate  unless  the  number  of  elements  in  the  y-direction  is  increased  to  a  point  where 
the  discontinuity  of  is  reduced  considerably.  This  would  result  in  a  very 

expensive  computational  analysis.  Therefore,  it  was  more  convenient  to  determine  cTj^ 
from  the  constitutive  relations. 
t.  Angle  —  ply  Laminate  [  ±  45], 

Since  the  stresses  were  determined  at  the  center  of  the  elements,  to  determine  the 
predicted  results  at  the  center  of  the  longitudinal  dimension  (x  =  L/2),  an  odd  number 
of  elements  were  used  when  using  Q8  or  Q9  element.  Because  these  elements  have 

midside  nodes,  displacements  were  also  directly  available  at  jt  =  .^.  In  the  case  of  04 

mr 

element,  two  analyses  had  to  be  carried  out.  To  obtain  the  solution  for  stresses,  an 
odd  number  of  elements  were  used.  However,  as  for  the  QS  and  ()9  elements,  to 
obtain  the  displacement  results  at  (L/2),  it  was  necessary  to  discretize  the  longitudinal 
direction  into  an  even  number  of  elements,  (odd  number  of  nodal  points)  so  that  a 
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set  of  nodal  points  would  be  located  at 


04  element 

Fig.  (12),  shows  plots  of  values  of  cr^  at  mid-surface  of  the  top  lamina  for  mesh 
refinement  in  the  y-direction.  The  number  of  elements  in  the  x-direction  was  kept 

constant  at  11  and  over  the  thickness,  each  layer  had  constant  rotation  about  the  x- 
and  y-axes.  As  observed  from  the  figure,  accuracy  of  the  axial  stress  is  unaffected  by 

refinement  in  the  \  direction.  However,  the  values  do  improve  somewhat  with 

refinement  in  the  longitudinal  (x)  direction  (Figs,  l,'^  and  14).  Fig.  (13)  shows  the 

efiect  of  refinement  along  x  into  21  and  31  elements  with  the  number  of  elements 

along  y  kept  constant  at  17.  Fig.  (14)  shows  the  effect  of  refinement  along  x  for  9 

elements  in  the  y-direction.  Fig.  (15)  shows  the  effect  of  thickness  refinement  on 
The  results  improve  significantly  with  refinement  near  the  free-edge.  The  same 
observation  can  be  made  for  in  Figs.  (16)  through  (18)  resjjectively.  The  calculated 
values  of  improve  slightly  with  refinement  in  the  x-direction  but  are  unaffected  by 
the  refinement  along  the  y  and  z-directions.  Accuracy  of  is  not  affected  by 
refinement  in  the  longitudinal  or  the  transverse  directions.  Figs.  (19)  through  (21). 

Accuracy  of  the  stresses  at  the  free-edge  improves  significantly  with  refinement 
over  the  thickness,  i.e.  as  N  is  increased  from  2  to  6.  This  refinement  does  not 

introduce  additional  elements  but  increases  the  degree  of  freedom  at  each  of  the  nodal 
points  by  assuming  rotation  to  be  constant  over  a  smaller  portion  of  the  thickness. 
Comparison  of  results  for  different  values  of  N  indicates  that  N=6,  in  comparison  to 
N=2  and  N=10,  gives  more  accurate  prediction  for  a.  Fig.  (15).  Further  refinement  to 
N=10,  in  comparison  to  N=2  and  N=6,  leads  to  better  estimation  of  Fig.  (22)  and 
T,.  Fig.  (23).  It  should  be  mentioned  that  thickness  refinement  does  not  influence  the 
accuracy  of  stres,ses  near  b=().  Stresses  at  the  free-edge,  however,  are  greatly 

influenced. 
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Figure  12:  X-stress  at  the  mid-surface  of  the  top  lamina  with  refinement 

in  y-direction;  Ang'e-ply  specimen. 
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Figure  13: 


X-stress  at  the  micl-surfnce  of  the  top  lamina  with  refinement 
in  >  d.rection:  An^le  plv  specimen. 


.20  2.40 


72 


+  Q4  ELEMENT  N=10  21X9 
A  04  ELEMENT  N=6  21X9 
O  04  ELEMENT  N=2  21X9 


°  "  PflGflNO  N=6 


Y/B 


Figure  15:  X-stress  at  midsurface  of  top  layer  with  thickness  refinement; 

Angle-ply  specimen. 
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Figure  16:  XY-stress  at  midsurface  of  top  layer  with  refinement 

y-direction;  Angle-ply  specimen. 
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Figure  17:  XY-stress  at  midsurface  of  top  layer  with  refinement  in 

x-direction:  A,  gle-ply  specimen. 
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Figure  20:  XZ-stress  at  z=h  with  refinement  in  ^-direction;  Angle-ply  spec¬ 

imen. 
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figure  23:  XZ-stress  at  z=h  with  thickness  refinement;  Angle-ply  specimen. 
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Comparison  of  accuracy  for  dif  fereni  types  of  elements 


Accuracy  of  the  stresses  calculated  using  different  types  of  elements  was  compared 
and  the  execution  time  for  each  case  noted.  Four,  eight,  and  nine-noded  elements  were 
used.  Two  different  mesh  sizes  viz.  9x9  and  21x9,  with  edge  elements  were  used 
with  four-noded  elements  (Fig.  6).  For  eight  and  nine-noded  elements  only  a  9x9 
mesh  with  'edge  elements',  i.e.  elements  along  the  edge  having  very  small  dimension  in 
the  y-direction,  was  used.  The  CPU  times  on  the  Cray  X-MP/28  for  this  set  of 
problems  are  shown  in  Table  (3).  Fig.  (24),  for  N=2,  shows  that  the  04  element  'vith 
the  21x9  mesh  predicts  the  axial  stress  at  the  free  edge  W’ith  11.47%  error  as  compared 
to  the  08  element  with  11.09%  and  09  element  with  10.99%.  For  N=t),  Fig.  (25).  th» 
error  for  the  04  element  reduced  to  2.68%,  for  08  to  2.22%,  and  for  09  to  2.12'? 

Fig.  (26),  for  N=10,  the  error  for  Q4  is  5.72%  and  for  08  is  5.11?o.  Throughout,  toc 
04  element  gives  significant  error  near  Y/b=0.  Comparison  of  inplane  shear  stress 
shows  that  the  results  at  the  free  edge  are  predicted  with  approximately  same  accuracy 
(Fig.  (27)  through  (29))  by  all  the  elements.  In  the  case  of  t„  (Fig.  (30)  through 
(32)),  Q4  with  9x9  and  21x9  mesh  refinement  gave  practically  the  same  results 
throughout  the  width  of  the  specimen  showing  that  the  system  is  rather  insensitive  to 
refinement  along  the  length  of  the  specimen. 

Further,  08  and  09  gave  the  same  result  at  the  free  edge  but  were  more  accurate 
than  04-  In  Fig.  (33),  the  longitudinal  displacement  predicted  at  X=L/2  and  the  top 
surface  of  the  specimen  using  08  element  is  shown  for  N=2,  N=6,  and  N=10.  In  Fig. 
(34)  and  (35)  corresponding  results  are  shown  for  N=2  and  N=6  in  the  case  of  Q9  and 
(J4  elements.  It  is  seen.  Fig.  (33),  that  for  N=2  the  displacements  are  slightly 
overpredicted  at  the  free  edge  and  for  N=6  and  N=10  displacements  are  slightly 
underpredicted.  Though  the  results  improve  with  refinement  from  N=2  to  N=6,  further 
improvement  with  refinement  to  N=l()  is  not  realized. 
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Tabu  3: 


Comparison  of  CPL  lime  on  X-MP/28  for  Q4,  Q8.  and  Q9  ele 
ments. 


N  Element  Mesh 

9X9 

04  - 

21X9 


CPU  on  X-MP/28 
(SEC) 

1.82 

4.479 


08 

9X9 

16.555 

09 

9X9 

28.181 

04 

9X9 

16.961 

21X9 

40.279 

08 

9X9 

190.222 

09 

9X9 

323.116 

04 

9X9 

60.755 

VjTT 

21X9 

145.421 

08 

9X9 

704.333 
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Figure  27:  XY-stress  for  N=2,  Q4,  Q8,  and  Q9  elements. 
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Figure  28:  XY-stress  for  N=6,  04,  Q8,  and  Q9  elements. 
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Figure  29:  XY-stresss  for  N-10,  Q4,  Q8,  and  g9  elements. 
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Figure  31:  XZ-stress  at  z=h  for  N=6,  Q4,  QS,  and  Q9  elements. 
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Figure  34:  Longitudinal  displacement  at  L’(L/2,y,2h)  of  09  element.  N=2, 

and  6. 
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Figure  35:  Longitudinal  displacement  at  L'(L/2.y.2h)  of  Q4  element,  N=2, 

and  6. 
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ti.  Cl  OSS- ply  laminale  [0 '90] 

Values  of  and  tr..  were  calculated  for  four-layer  crass-ply  [0/90\  using 
equations  (2.33)  and  (2..34)  as  was  done  in  the  case  of  t,.  for  the  angle-ply  laminate 
[±45], 

.As  shown  in  Fig.  (36)  mesh  refinement  in  the  y-direction  does  not  lead  to  any 
noticeable  improvement  in  the  calculated  values  of  Fig.  (37)  indicates  that 

appro.\imation  of  is  not  dependent  on  axial  refinement  either.  Results  obtained  for 
refinement  through  the  thickness.  Fig.  (38),  for  11x17  mesh,  17  being  the  number  of 
elements  along  the  y-direction,  indicate  tiiat  tin:,  refinement  profoundly  influences  the 
results.  However,  the  traction-free  biiundary  condition  at  the  free  edge  is  not  satisfied 
and,  therefore,  the  solution  does  not  match  Pagano's  results.  The  calculated  valt 
cr^.  are  in  error  up  to  3()%  even  with  refinement,  for  y^b  equal  to  0.8.  Near  the 
free-edge,  the  error  is  quite  large.  The  values  of  ajj  determined  from  equation  (37) 
were  either  exactly  zero  or  close  to  zero.  This  is  clearly  wrong  and  repre.sents  a  setiOus 
limitation  of  the  theory. 


4.1.4  Analysis  of  22-Layer  Free-Edge  Delamination  Specimen 

The  procedure  developed  was  applied  to  a  22-layer  coupon  with  fiber  orientation 
of  [( 25.5 /  — 25.5  )j /  90],  which  was  previously  solved  by  Chang  [ll]  and  Dandan  [l2]. 
The  laminate  width  and  ply  thickness  were  taken  as  1.0  inch  and  0.00505  inch, 
respectively.  The  material  properties  of  a  lamina  were  the  same  as  used  in  previous 
investigations  [ll],  [12]. 


=  19.26*10*’  ipsi) 

^22  =  -^33  =  1-32*  lO**  (psi) 
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YZ-stress  at  z=h  with  mesh  refinement  in  x-direction;  Cross-ply 
specimen. 


Due  to  the  symmetry  oi  the  laminate,  it  was  only  neces-sary  to  consider  11  lamina  in 
the  analysis.  Bitsed  on  the  analysis  performed  on  the  four-layer  coupon  specimen,  it 
was  concluded  that  refinement  along  the  x-direction  improves  the  accuracy  of  the 

results  more  effectively  than  the  refinement  along  the  y-direction.  For  this  reason,  the 
22-layer  specimen  was  discretized  into  22  elements  in  the  x-direction  and  13  in  the 
y-direction.  All  elements  have  the  same  dimensions  in  the  x-direction.  However,  in 
the  y-direction,  the  edge  elements  have  a  width  of  ().(X)5  inch,  with  the  remaining  por¬ 
tion  discretized  shown  in  IMg.  (6).  The  formulations  in  [ll]  and  [12]  w'ere  based  on  a 
two-dimensional  nuxiel  dependent  on  y  and  /  cixirdinates  only.  In  [12],  22  elements 
were  used  in  the  /.-direction  and  14  in  the  y-direction.  In  the  present  investigation, 
the  discreii/Jition  is  along  the  x-  and  y-coordinates.  Therefore,  it  was  not  feasib  to 
match  the  mesh  with  that  used  in  [ll]  and  [12]. 

To  compare  the  results  with  those  given  in  [12]  and  [ll],  two  different  analyses 
were  performed.  In  [12]  the  results  are  given  at  the  center  of  each  layers.  In  the 

present  investigation,  components  are  calculated  at  the  centers  of  the  laminae,  and 

based  upon  the  effectiveness  of  the  refinement  over  the  thickness  noticed  in  the 

analysis  of  the  four-layer  specimen,  each  layer  is  divided  into  three  sublayers,  forming 
a  total  of  33  layers.  Cj,  is  calculated  at  the  interface  in  the  present  investigation. 
However,  in  order  that  these  results  be  comparable  to  those  obtained  in  [12]  and  [ll] 
at  the  center  of  the  lamina,  each  layer  must  be  subdivided  into  two  sublayers,  forming 
a  total  of  22  layers,  producing  the  necessary  interfaces.  An  average  axial  strain  of 
(0.95414* lO"*)  was  applied,  which  is  the  same  as  applied  by  Chang  [ll]. 

Fig.  (39)  shows  the  cross-section  of  the  upper  half  of  the  symmetric  laminate  and 
defines  the  location  for  plotting  of  results.  Fig.  (40)  and  Fig.  (41)  show  the 
distribution  of  <x^  along  the  center  of  the  11th  and  the  5th  layer  for  different  (y/B) 
ratios.  The  results  obtained  from  I'SDT  match  those  obtained  from  the  higher  order 


HX) 


element  [ll]  and  the  axisymmetric  model  [12].  Similarly,  Fig.  (42)  shows  that  c,  at 

the  center  of  the  11th  layer  can  be  obtained  accurately,  though  the  concentration  of 

the  stress  at  the  free-edge  does  not  match  that  obtained  by  the  higher  order  element 
[ll]  Distribution  of  along  R^^,  as  shown  in  Fig.  (43),  closely  follows  the  result 

given  in  [ill  However,  the  stresses  determined  from  the  FSDT  are  in  considerable 
error  near  the  free  edge.  Figs.  (44)  to  (47)  show  the  plots  for  cr^j  along  Rl,  R5,  R6, 
Rll,  and  Figs.  (48)  to  (50)  indicate  the  cr,.  along  R5.  R6,  and  Rll.  It  is  evident  that 
both  cr^.  and  <7^.  are  predicted  reasonably  well  by  the  FSDT  in  the  regions  y/B<().l. 
However,  for  y  B>0.9,  the  predicted  stresses  for  O’,,  even  result  in  signs  different 

from  those  in  [ill  Results  for  O’,,,  do  not  satisfy  the  free-edge  stress  at  y/B  =  1.0. 
Figures  (51)  through  (53)  show  the  through-the-thickness  distribution  of  o’..,  cr,.,  and 
O’,,  along  the  free-edge.  The  results  from  the  FSDT  do  not  agree  with  those  from  [ll] 
and  [12]  and  are,  apparently,  quite  wrong.  Furthermore,  Figs.  (54)  to  (57)  show  the 
distribution  of  <r..  along  Rl,  R5,  R6,  and  Rll.  The  stress  distribution  obtained  from 
the  present  approach  is  quite  different  from  that  given  by  [l2l 
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Figure  47’.  ZX -stress  at  Rll  for  22-layer  coupon. 
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Figure  53:  Through  the  thickness  distribution  of  YZ-stress  at  y/b-0.995  for  the 

22-layer  coupon. 
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FREE  EDGE 
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FREE  EDGE 


4.1^  Summary  and  Conclusions 

1.  The  Q8  and  Q9  elements  are  superior  to  the  Q4  element  for  the  analysis  of 

plates.  However,  for  analysis  of  coupons  under  uniform  extension,  the  per¬ 
formance  of  the  Q8  and  Q9  elements  was  comparable  to  the  Q4  element. 
The  CPU  times  needed  for  the  Q8  and  Q9  were  much  greater  than  required 

for  the  analysis  using  Q4  elements.  Thus,  the  Q4  element  appears  to  give  the 

better  combination  of  accuracy  and  economy  of  computational  effort  for  cou¬ 
pons. 

2.  To  improve  the  predicted  results  for  cr^  and  t,,  at  y/b=0,  refinement  along 

the  longitudinal  direction  seems  to  be  most  effective.  However,  the  solution 

at  the  free-edge  depended  upon  thickness  refinement.  The  best  results  at  the 

free-edge  for  <r^  were  for  N=6,  while  for  they  were  for  N=10. 

3.  The  results  for  r^.  and  t^.  did  not  improve  significantly  with  refinement  in 
X-  or  y-directions.  The  thickness  refinement  improved  the  predicted  results  at 
the  free-edge  for  t^.,  but  did  not  satisfy  the  traction-free  boundary  condition 
at  the  free-edge  for  t,.. 

4.  The  predicted  results  for  cr,  are  all  either  zero  or  close  to  zero  for  the  four- 
layer  delamination  specimens.  However,  for  the  22-layer  delamination  speci¬ 
men,  they  were  oscillatory  along  the  z-coordinate  and  quite  different  from  the 
results  obtained  by  Chang  [ll]  and  Dandan  [12]. 

5.  The  predicted  displacements  for  free-edge  delamination  specimens  are  close  to 
Pagano's  results.  For  N=2  the  displacement  results  are  slightly  overpredicted, 
and  for  N=6  and  N=10,  the  displacements  are  slightly  underpredicted. 
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Section  V 


DISCUSSION 

In  the  present  research  program,  the  theoretical  studies  carried  out  have  included  \ 

two  distinct  approaches  to  the  problem  of  stress  analysis  of  composite  laminates.  One 
approach  consisted  of  a  specialization  of  the  three-dimensional  elastostatics  theory  to  the 
case  of  free-edge  delamination  specimens  under  uniform  axial  strain  in  which  the  stress 
field  is  independent  of  the  longitudinal  coordinate.  The  other  approach  consi.sted  of 
development  and  application  of  theories  of  laminated  plates  to  the  problem  of  fre.  euge 
delamination.  The  present  report  covers  study  of  the  applicability  of  the  existing 
’’discrete  laminate  theory”  as  one  of  the  first  steps  in  the  research  program.  This 
investigation  has  served  to  establish  the  pattern  for  later  efforts  on  improved  theories 
for  study  of  free-edge  delamination  of  composite  laminates. 

The  equations  of  the  discrete  laminate  theory  are  well  known.  That  have  been 
included  in  this  report  for  completeness  and  for  the  purpose  of  pnintino  out  the 
theoretical  assumptions  inherent  in  the  theory.  Mawenya  and  Davis  [S]  had  previously 
presented  a  finite  element  implementation  of  the  theory,  but  they  provided  no  details 
and  did  not  apply  the  theory  to  free-edge  delamination. 

The  equations  of  the  discrete  laminate  theory,  essentially  treating  a  laminate  plate 
as  a  stacking  of  Mindlin  plates,  assumes  linear  variation  of  the  in-plane  displacements 
over  the  thickness  of  each  layer  ensuring  continuity  of  displacements  at  the  interfaces. 

As  part  of  the  present  research  effort,  the  equations  have  been  written  in  matrix  form 
such  that  the  matrix  of  operators  is  self-adjoint  in  an  appropriate  linear  vector  qiace. 

The  general  formulation  was  written  for  the  dynamic  problem  in  the  convolution 
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product  space  and  then,  for  the  elastostatic  problem,  a  special  form  was  defined  in  the 
inner  product  space.  For  such  systems  of  equations,  standard  techniques  are  available 
for  the  construction  of  variational  principles  and  for  identification  of  consistent 
boundary  operators.  General  variational  formulations  with  extensions  and  useful 
specializations  have  been  explicitly  developed  for  the  problem.  One  specialized 
variational  formulation,  corresponding  to  the  popular  potential  energy  theory,  has  been 
used  to  develop  finite  element  procedures. 

Three  different  isoparametric  finite  element  interpolation  schemes  viz.,  the 
four-point  Lagrangian,  the  nine-point  Ijgrangian,  and  the  eight-point  serendipity,  have 
been  implemented  in  a  computer  program  written  initially  for  an  IBM  3081  mainframe 
computer  and  later  modified  to  run  on  a  CRAY-XMP/28.  These  finite  element 
procedures  were  verified  through  application  to  homogeneous  as  well  as  sandwich  plate 
problems  for  which  solutions  are  available.  Their  effectiveness  in  modelling  the 
stress-distribution  in  free-edge  delamination  specimens  has  been  examined. 

The  present  investigation  indicated  that  the  discrete  laminate  theory  of  laminated 
plates  is  quite  effective  in  modelling  displacements  in  plates  subjected  to  arbitrary 
transverse  loads.  The  shear  effects  can  be  allowed  for  satisfactorily.  For  such 
problems,  the  higher  order  elements  viz.,  the  nine-point  Lagrangian  and  the  eight-point 
serendipity,  performed  better  than  the  simple  four-point  Lagrangian  element  However, 
for  application  to  free-edge  delamination,  the  entire  approach  is  inadequate.  The  stress 
distribution  obtained  for  the  example  problems  was  reasonably  good  with  respect  to 
in-plane  stresses,  but  the  theory  could  not  give  reasonable  estimates  of  the  other  three 
components  of  stress.  The  traction-free  edge  condition  could  not  be  modelled.  The 
stresses  had  to  be  calculated  directly  from  the  stress-strain  relationships  because  use  of 
equilibrium  equations  for  determination  of  shear  stress  and  transverse  stress  would 
involve  numerical  differentiation  of  quantities  for  which  estimates  only  at  a  finite 


number  of  points  were  available.  Mesh  refinement  to  get  a  sufficiently  large  number 
of  points  would  make  the  cost  of  analysis  prohibitive. 

The  studies  showed  that  refinement  along  the  length  or  the  width  of  the  specimen 
had  relatively  little  effect  on  the  quality  of  results.  Refinement  over  the  thickness 
i.e.,  subdivision  of  each  lamina  into  a  number  of  sublayers,  helped  improve  accuracy. 
This  suggests  that  a  distribution  of  in-plane  stresses  of  an  order  higher  than  linear 
over  the  thickness  of  each  lamina  might  represent  the  actual  stress  distribution  more 
closely. 

In  order  to  satisfy  the  traction-free  conditions  along  the  free-edges,  it  is  necessary 
that  edge  tractions  appear  as  field  variables  in  the  set  of  field  equations.  .4n 
alternative,  of  course,  is  to  use  Lagrange  multiplier  techniques  to  enforce  constraint 

For  direct  use  of  equilibrium  equations  to  determine  the  shear  stresses  and  the 
direct  transverse  stress,  explicit  introduction  of  interface  tractions  as  field  variables  in 
the  theory  would  avoid  the  need  for  expensive  numerical  differentiations.  Also  this 
would  ensure  continuity  of  traction  across  interfaces  and  perhaps  yield  better 
approximation  for  the  interfacial  stresses. 

The  discrete  laminate  theory  discussed  in  this  report  is  based  on  assumptions 
regarding  transverse  and  in-plane  displacements.  An  alternative  is  to  assume  variation 
of  in-plane  stresses  and  to  derive  the  other  stress  components  through  equilibrium 
equations.  If  force  resultants  appear  in  the  expressions  for  stresses,  and  are  regarded  as 
field  variables  of  the  problem,  the  constitutive  relationships  for  these  need  to  be 
established.  If  it  is  assumed  that  there  is  no  interfacial  slip,  it  would  be  impossible 
for  any  layer  to  deform  independently  of  the  others.  This  would  necessarily  lead  to  a 
coupling  in  the  constitutive  relations  for  the  force  resultants  of  individual  layers. 

Pagano’s  [lO]  theory  which  uses  the  assumption  of  linear  variation  of  in-plane 
stress  components  over  the  thickness  of  each  layer  or  sublayer,  satisfies  equilibrium 
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pointwise,  satisfies  constitutive  equations  and  inienacial  continuity  of  tractions  as  well 
as  displacements,  and  can  satisfy  uaction  boundary  conditions  for  free-edge  delamination 
specimens  exactly,  appears  to  be  an  appropriate  approach  for  determination  of  stress 
fields  in  composite  laminate  plates  with  natural  boundary  conditions.  The  case  of 
free-edge  delamination  specimens  is  a  specialization  of  the  general  theory.  The  theory 
has  been  difficult  to  use  because  of  the  large  number  of  field  variables  involved  and 
limitations  on  computational  capabilities.  A  research  effort  directed  towards 
development  of  finite  element  models  based  on  Pagano's  theory  or  development  of  other 
methods  of  solution  of  the  set  of  differential  equations  could  be  useful. 
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Appendix  A 

VARIATIONAL  FORMULATION 


Often,  obtaining  an  approximate  solution  to  a  coupled  boundary  value  problem 
relies  on  appropriate  variational  formulation.  Following  Sandhu's  [6]  [7]  [8]  [9], 
extension  of  Mikhlin's  [l3]  basic  variational  theorem  to  coupled  linear  boundary  value 
problems  including  nonhomogenous  boundary  condition,  we  present  here  a  summary  of 
the  basic  concepts  for  setting  up  the  variational  formulation  applicable  to  the  problem 
of  laminated  plates. 


A.1  PRELIMINARIES 


A.1.1  Boundary  Value  Problem 

Consider  the  boundary  value  problem 

Au  =  /  on  jR  (A.l) 

Cu  =  g  on  QR  (A.2) 

where  is  the  boundary  of  the  open  connected  region  R  in  an  euclidean  space.  R  is 
the  closure  of  R.  A  and  C  are  linear  bounded  operators.  Let  and  be  linear 
vector  spaces  defined  on  the  regions  indicated  by  the  subscripts,  and  Wj,,W^  be  dense 
subsets  in  Vj,  and  respectively.  Then  the  differential  operators  A  and  C  can  be 
regarded  as  the  transformations 

AiW^^Vj^  (A.3) 


C:W 


9f! 
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A.1.2  Bilinear  Mapping 

Let  V  and  S  be  linear  vector  spaces.  A  bilinear  mapping  B:VxV-*S  assigns  to 
each  ordered  pair  of  vectors  m,v€V  an  element  in  S.  Furthermore,  bilinearity  is 
satisfied  for  u,,Uj,v,,v^,u,v  €  V,  if 


■V 

Kou,  +M2.V)  =  aB(.u^,v)  +  Biu^yV) 

(A.4) 

; 

B(u,0(V|  +  V2)  =  aB(u,v^)  + Biu.v^) 

(A.5) 

where  a  is  scalar.  For  convenience,  we  shall  use  the  notation. 

Bg{u,v)  =  <  u.v  > ^  (A.6) 

To  set  up  a  variational  formulation,  symmetric,  nondegenerate  bilinear  mappings  are 
used,  i.e., 

<u,v>  <v,u>^  (A.7) 

and 

<«,v>  =  0  for  all  V  if  and  only  if  u-i)  (A.8) 

A.U  Self-Adjoint  Operator 

An  operator  A*  on  V  is  said  to  be  the  adjoint  of  A  with  respect  to  symmetric 
bilinear  mapping  Bg:VxV-*S,  where  S  is  a  linear  vector  space,  if 

<  «,  Av  =  <  V,  A  u>ju  + Z)^(  v,m)  (A.9) 

for  all  u  and  v  €  V  and  where  represents  quantities  associated  with 

boundary  of  R.  If  A=  A*,  then  A  is  said  to  be  self-adjoint.  If  A  is  a  self-adjoint 

operator,  then  D^iv,u)  is  antisymmetric,  i.e., 

D^(v,u)  = -/7^jp(u,v)  (A.IO) 

Furthermore,  A  is  said  to  be  symmetric  wdth  respect  to  the  bilinear  mapping,  if 

<u,Av>j^=  <v,Au>g  (A.ll) 

The  boundary  operator  ('  is  said  to  be  consistent  with  the  self-adjoint  operator  A  if 
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Z)g^(v.u)  =  <  u.Cv  -  <  v,Cu 


(A.12) 


A.1.4  Gateaux  Differential 

If  Cl:V-*S.  where  V  is  sucli  that  if  u,  u6\  .  ut  AuCV’  for  scalar  X, 
the  Gateaux  differential  o(  fl  (u)  along  a  path  u  is  defined  by 

S_  n(u)  =  (A.13) 

"  x-o  A 

where  u  is  referred  tn  as  the  path. 

A.2  THEOREM 

Tor  the  field  equations  (A.l )  we  define 

Cl(u)  =  <u,Au>  ^,  —  2<u  ,f>  (A.M) 

The  Gateaux  differential  of  fl  is: 

8  fl(  )  <  u  +  \u,A(u  -f  Xu)>  —2<u  +  \u,f>  —  <u,Au>+2<  u,f> 

=  <u,AU>  A  <U,Au>  —2  <u,f> 

-  2  <Q,Au  —  f> 

The  Gateaux  differential  vanishes  at  the  solution  «  =  u,.  where  Au^  —  f  =  0.  Conversly, 
if  SgflCu)  vanishes  for  all  U,  nondegeneracy  of  <  ,  >  implies  Au„  —  f  =  0  If  the  range 
of  the  bilinear  mapping  is  the  real  line,  vanishing  of  the  function  O  would  imply  its 
minimum,  maximum,  or  stationary  value,  depending  upon  the  operator  A  being  positive, 
negative  or  semi-definite. 


* 


1-^0 


A3  UNEAR  COUPLED  PROBLEMS 


% 


The  above  discussion  for  a  single-valued  function  u  can  be  extended  to  the  case  of 
several  variables.  If  there  are  n  variables,  V  is  defined  as  the  direct  sum 

V  =  V,+V,  + - +  V,  (A.16) 

and  an  element  u€V  is  an  n-tuple  (u,  .Uj, _ ,u„)  with  u,  €  for  i=l,  2,  n.  A 

bilinear  mapping  on  V  is  defined  as 

<u,v>  =  <u  ,v  >  +<u,,v,>  + . +  <u  ,v>  _  (A.17) 

where  <  ,  is  defined  for  components  of  respectively. 

If  the  field  and  boundary  condition  of  a  linear  coupled  boundary  value  problem 
are: 

ft 

£a,  .u^  =  /.  on  R  (A.18) 

y-l 

n 

f-l,2,..Vt  (A.19) 

J-i 

the  governing  functional  based  on  Hqs.  (A.18)  and  (A.19)  is 

n  ti  u  n 

0  («)  =  £  <  u, ,  J:  A,  >  ,  +  £  <  “,  •  £  C,^  -  2s,  >  „  (AJO) 

I- 1  j~  1  »- 1  y- 1 

The  set  of  operators  A,^  is  said  to  be  self-adjoint  with  respect  to  the  bilinear  mapping, 
if 


<  V.,  7*  A  M  >  -~  <  7*  u  ,  A  v  >  „  +  D-.(u  ,v) 

tj  j  f  /  ji  I  Jf  dP  /'  I 


(A.21) 


J-i 


where  D^Uj,v)  represents  quantities  associated  with  boundary  QR  of  R.  The  boundary 
operators  are  said  to  be  consistent  with  the  field  operator  if 


D.lu  ,v)  -  ^<u  ,C  V  >  <  V  ,y'  C.  u  > 

j  *  j  p  *  V  J 


(A32) 


I 


;  1 


1.^1 


Appendix  B 


SOLUTION  OF  SANDWICH  PLATE 


B.1  PRELIMINARIES 

Based  on  I’lanlema  [14],  ilie  series  solution  for  sandwich  plate  was  calculated.  In 
this  case,  the  lx.‘ndinj;  of  plates  is  assumed  to  l)c  due  to  bending  of  stiff  layers  and  the 
shear  deformation  of  the  core  layer,  so  that 

w  =  +  w  (B.1 ) 

h  i 

where  the  transverse  load  is  the  only  applied  load  and  the  stiff  layers  are  isotropic,  vv^ 
satisfies  the  equation; 

=  q  (a2) 


where 


Et\ 

D  =  ' 


12(1-1/") 
and  w,  satisfies  the  equation: 


where 


-SV  w  =  q 


.Kill 


(B.3) 


1.12 


B.2  METHOD  OF  SOLUTION 

For  a  simply  supported  rectangular  plate  the  transverse  displacement  due  to 
bending  may  be  represented  by  Fourier  series  in  the  form: 

w(x,y)  =  £  2^  a^sin 

m-l  n-1 

and  corresponding  load  by: 


«?(>•■.>)  =  Z  Z ^  ^ 


Multiplying  both  sides  of  (B.5)  by  sin  sin  -  and  integrating  over  the  domain, 

a  0 

for  q(x,y)=  q„. 


^mn  2 


Substituting  (B.4)  and  (R5)  into  (B,2)  and  evaluating  a„,  w^ix,y)  is: 

„ Xxj.)  =  £ f 

rr^D  „  „  +  (n/&)^]^ 

Similarly,  (E3)  can  be  solved  to  yield: 


.  16^0  y«  Y»  sin(m7rj:/a)sin(niry/&) 

*  v*S  „  „  mnl(m/a)^  +  (n/W^] 


(m,n  =  1,3,...) 


The  series  solution  was  obtained  through  summation  of  151  terms  in  the  case  of 
isotropic  sandwich  plate  for  a>>=b>10  inches  and  =  1.0  lb/in\  In  the  orthotropic  case 
the  results  given  by  Mawenya  and  Davis  [5]  were  used  for  comparison  purposes. 
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